Show code cell content Hide code cell content

# Source the package setup script
source("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/scripts/00_setup_packages.R");

The C++ toolchain required for CmdStan is setup properly!

* Installing CmdStan from https://github.com/stan-dev/cmdstan/releases/download/v2.32.1/cmdstan-2.32.1.tar.gz

Warning message:
"An installation already exists at C:\Users\bmc82/.cmdstan/cmdstan-2.32.1. Please remove or rename the installation folder or set overwrite=TRUE."

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

Loading required package: usethis

Linking to libssh v0.11.0

Attaching package: 'kableExtra'

The following object is masked from 'package:dplyr':

    group_rows

Warning message:
"package 'FSA' was built under R version 4.4.3"

Registered S3 methods overwritten by 'FSA':
  method       from
  confint.boot car 
  hist.boot    car 

## FSA v0.9.6. See citation('FSA') if used in publication.
## Run fishR() for related website and fishR('IFAR') for related book.

Attaching package: 'plotly'

The following object is masked from 'package:ggplot2':

    last_plot

The following object is masked from 'package:stats':

    filter

The following object is masked from 'package:graphics':

    layout

Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

Loading required package: rJava

Loading required package: leaps

Loading required package: sp

Attaching package: 'raster'

The following object is masked from 'package:plotly':

    select

The following object is masked from 'package:dplyr':

    select

Attaching package: 'recolorize'

The following object is masked from 'package:plotly':

    add_image

Warning message:
"package 'webshot2' was built under R version 4.4.3"

This is bayesplot version 1.11.1

- Online documentation and vignettes at mc-stan.org/bayesplot

- bayesplot theme set to bayesplot::theme_default()

   * Does _not_ affect other ggplot2 plots

   * See ?bayesplot_theme_set for details on theme setting

Attaching package: 'bayesplot'

The following object is masked from 'package:plotrix':

    plot_bg

This is posterior version 1.6.0

Attaching package: 'posterior'

The following object is masked from 'package:bayesplot':

    rhat

The following objects are masked from 'package:raster':

    %in%, match

The following objects are masked from 'package:stats':

    mad, sd, var

The following objects are masked from 'package:base':

    %in%, match

Attaching package: 'gridExtra'

The following object is masked from 'package:dplyr':

    combine

Loading required package: Rcpp

Loading 'brms' package (version 2.22.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').

Attaching package: 'brms'

The following object is masked from 'package:bayesplot':

    rhat

The following object is masked from 'package:stats':

    ar

Attaching package: 'scales'

The following object is masked from 'package:plotrix':

    rescale

── Attaching core tidyverse packages ──── tidyverse 2.0.0 ──
✔ forcats   1.0.0     ✔ readr     2.1.5
✔ lubridate 1.9.4     ✔ stringr   1.5.1
✔ purrr     1.0.2     

── Conflicts ────────────────────── tidyverse_conflicts() ──
✖ readr::col_factor()      masks scales::col_factor()
✖ gridExtra::combine()     masks dplyr::combine()
✖ purrr::discard()         masks scales::discard()
✖ raster::extract()        masks tidyr::extract()
✖ plotly::filter()         masks dplyr::filter(), stats::filter()
✖ kableExtra::group_rows() masks dplyr::group_rows()
✖ dplyr::lag()             masks stats::lag()
✖ raster::select()         masks plotly::select(), dplyr::select()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Attaching package: 'patchwork'

The following object is masked from 'package:raster':

    area

Attaching package: 'cowplot'

The following object is masked from 'package:patchwork':

    align_plots

The following object is masked from 'package:lubridate':

    stamp

The following object is masked from 'package:ggpubr':

    get_legend

This is loo version 2.8.0

- Online documentation and vignettes at mc-stan.org/loo

- As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session. 

- Windows 10 users: loo may be very slow if 'mc.cores' is set in your .Rprofile file (see https://github.com/stan-dev/loo/issues/94).

Attaching package: 'extraDistr'

The following object is masked from 'package:purrr':

    rdunif

The following objects are masked from 'package:brms':

    ddirichlet, dfrechet, pfrechet, qfrechet, rdirichlet, rfrechet

Warning message:
"package 'leaflet' was built under R version 4.4.3"

Loading required package: StanHeaders

rstan version 2.32.6 (Stan version 2.32.2)

For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
change `threads_per_chain` option:
rstan_options(threads_per_chain = 1)

Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file

Attaching package: 'rstan'

The following objects are masked from 'package:posterior':

    ess_bulk, ess_tail

The following object is masked from 'package:raster':

    extract

The following object is masked from 'package:tidyr':

    extract

Warning message:
"package 'ggdist' was built under R version 4.4.3"

Attaching package: 'ggdist'

The following objects are masked from 'package:brms':

    dstudent_t, pstudent_t, qstudent_t, rstudent_t

Attaching package: 'ape'

The following objects are masked from 'package:raster':

    rotate, zoom

The following object is masked from 'package:glmulti':

    consensus

The following object is masked from 'package:ggpubr':

    rotate

The following object is masked from 'package:dplyr':

    where

Bioconductor version '3.20' is out-of-date; the current release version '3.21'
  is available with R version '4.5'; see https://bioconductor.org/install

Attaching package: 'BiocManager'

The following object is masked from 'package:devtools':

    install

treeio v1.30.0 Learn more at https://yulab-smu.top/contribution-tree-data/

Please cite:

LG Wang, TTY Lam, S Xu, Z Dai, L Zhou, T Feng, P Guo, CW Dunn, BR
Jones, T Bradley, H Zhu, Y Guan, Y Jiang, G Yu. treeio: an R package
for phylogenetic tree input and output with richly annotated and
associated data. Molecular Biology and Evolution. 2020, 37(2):599-603.
doi: 10.1093/molbev/msz240

Attaching package: 'treeio'

The following object is masked from 'package:raster':

    mask

Attaching package: 'TreeTools'

The following object is masked from 'package:treeio':

    MRCA

ggtree v3.14.0 Learn more at https://yulab-smu.top/contribution-tree-data/

Please cite:

Guangchuang Yu. Using ggtree to visualize data on tree-like structures.
Current Protocols in Bioinformatics. 2020, 69:e96. doi:10.1002/cpbi.96

Attaching package: 'ggtree'

The following object is masked from 'package:TreeTools':

    MRCA

The following object is masked from 'package:ape':

    rotate

The following objects are masked from 'package:raster':

    flip, rotate

The following object is masked from 'package:tidyr':

    expand

The following object is masked from 'package:ggpubr':

    rotate

Cluster Analysis

---

Question

Do camouflage males and females exhibit discrete groupings in dorsal body color pattern variables that suggest distinct morph types?

---

Objective

Test for clustering of dorsal body color pattern metrics (e_max, Filter_max, e_prop, R, G, B).

---

Method

1. Load cleaned data.

We start by loading the cleaned data from the “01_data_cleaning” pipeline. This data has already undergone transformations and contains relevant metrics for our models.

data_c1_clean <- read.csv("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/cleaned/data_c1_clean.csv")

data_c1_clean <- data_c1_clean %>%
  mutate(BCPD = paste0("BCPD_", sprintf("%04d", as.numeric(BCPD))))

---

2. Prepare data for clustering.

Individuals with distinct color patterns are visually classified and removed prior to cluster analyses. Also, male and female data are separated to eliminate Sex as a confounding variable.

# First, subset male and female data
data_c1_F = data_c1_clean[data_c1_clean$Sex == 'F',]
data_c1_M = data_c1_clean[data_c1_clean$Sex == 'M',]


# Subset color types to cluster
# First, make dataframe out of known morphs so we can append these back later. 
data_c1_F_morphs <- subset(data_c1_F, Morph_original %in% c("fluor","beige_flour","mimic_purple","mimic_orange_type1",
      "mimic_orange_type2","green","xmas","brown_beige_border"))

data_c1_M_morphs <- subset(data_c1_M, Morph_original %in% c("fluor","beige_flour","mimic_purple","mimic_orange_type1",
      "mimic_orange_type2","green","xmas","brown_beige_border"))


# Then, remove these morphs from the dataset that will be clustered
data_c1_F <- data_c1_F %>%
    filter(!(Morph_original %in% c(
      "fluor","beige_flour","mimic_purple","mimic_orange_type1",
      "mimic_orange_type2","green","xmas","brown_beige_border")))

data_c1_M <- data_c1_M %>%
    filter(!(Morph_original %in% c(
      "fluor","beige_flour","mimic_purple","mimic_orange_type1",
      "mimic_orange_type2","green","xmas","brown_beige_border")))

Categorical variables

For categorical and binary predictors, R automatically dummy-codes the variables when they are stipulated as factors. Thus, we don’t need to worry about scaling them. Also, for our purposes, we are not using any categorical variables for clustering. We will just convert the categorical variables in the dataframe to factors to be consistent.

Continuous variables

We standardize continuous predictors by centering and scaling them (dividing by two standard deviations). Standardizing variables before clustering helps ensure that all features contribute equally to the clustering process, leading to more meaningful and unbiased clusters. It is a fundamental preprocessing step, especially when variables have different scales or units.

# Convert categorical variables to factors
columns_to_convert_c1 <- c("Sex", "Morph_original")

data_c1_F <- data_c1_F %>%
    mutate(across(all_of(columns_to_convert_c1), as.factor))

data_c1_M <- data_c1_M %>%
    mutate(across(all_of(columns_to_convert_c1), as.factor))


# Convert continuous variables to numeric
columns_to_convert_c1 <- c("e_max", "Filter_max", "e_prop", "R_c", "G_c", "B_c")

data_c1_F <- data_c1_F %>%
    mutate(across(all_of(columns_to_convert_c1), as.numeric))

data_c1_M <- data_c1_M %>%
    mutate(across(all_of(columns_to_convert_c1), as.numeric))

# Standardize continuous variables
columns_to_scale <- c("e_max", "Filter_max", "e_prop", "R_c", "G_c", "B_c")

# Scale the specified columns
scaled_columns_F <- as.data.frame(scale(data_c1_F[, columns_to_scale]), scale = 2 * sd(data_c1_F))

scaled_columns_M <- as.data.frame(scale(data_c1_M[, columns_to_scale]), scale = 2 * sd(data_c1_M))

# Combine scaled columns with unscaled columns
data_c1_F_scaled <- cbind(
  data_c1_F[, !(colnames(data_c1_F) %in% columns_to_scale)], #selects all columns from data_c1_clean except the ones specified in columns_to_scale.
  scaled_columns_F
)

data_c1_M_scaled <- cbind(
  data_c1_M[, !(colnames(data_c1_M) %in% columns_to_scale)], #selects all columns from data_c1_clean except the ones specified in columns_to_scale.
  scaled_columns_M
)

---

3. Principal Components Analysis (PCA)

PCA is conducted on scaled data to reduce dimensionality and identify key components.

data_c1_F_scaled.pca <- prcomp(data_c1_F_scaled[, 5:10], scale. = TRUE)
data_c1_M_scaled.pca <- prcomp(data_c1_M_scaled[, 5:10], scale. = TRUE)

Show code cell source Hide code cell source

# Adjust scree plot for females with custom y-axis label
PCA_screeplot_Female <- fviz_screeplot(data_c1_F_scaled.pca, addlabels = TRUE, ggtheme = theme_bw(base_size = 10)) +
  ggtitle("Scree Plot") +
  labs(y = "Explained Variance (%)") +  # Custom y-axis label
  scale_y_continuous(breaks = seq(0, 100, by = 10), labels = paste0(seq(0, 100, by = 10), "%")) +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.title.y = element_text(size = 12, margin = margin(r = 10)),
    axis.text.y = element_text(size = 10)
  )

PCA_screeplot_Male <- fviz_screeplot(data_c1_M_scaled.pca, addlabels = TRUE, ggtheme = theme_bw(base_size = 10)) +
  ggtitle("Scree Plot") +
  labs(y = "Explained Variance (%)") +  # Custom y-axis label
  scale_y_continuous(breaks = seq(0, 100, by = 10), labels = paste0(seq(0, 100, by = 10), "%")) +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.title.y = element_text(size = 12, margin = margin(r = 10)),
    axis.text.y = element_text(size = 10)
  )

# Adjust variable contribution plot for PC1 with custom y-axis label
fviz_cos2_F_PC1 <- fviz_cos2(data_c1_F_scaled.pca, choice = "var", axes = 1, ggtheme = theme_bw(base_size = 10)) +
  ggtitle("Variable Contributions to PC1") +
  labs(y = "Contribution (%)") +  # Custom y-axis label
  scale_y_continuous(breaks = seq(0, 1, by = 0.1), labels = paste0(seq(0, 100, by = 10), "%")) +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.title.y = element_text(size = 12, margin = margin(r = 10)),
    axis.text.y = element_text(size = 10),
    plot.margin = margin(t=10, b=100, l = 10, r = 10)
  )

fviz_cos2_M_PC1 <- fviz_cos2(data_c1_M_scaled.pca, choice = "var", axes = 1, ggtheme = theme_bw(base_size = 10)) +
  ggtitle("Variable Contributions to PC1") +
  labs(y = "Contribution (%)") +  # Custom y-axis label
  scale_y_continuous(breaks = seq(0, 1, by = 0.1), labels = paste0(seq(0, 100, by = 10), "%")) +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.title.y = element_text(size = 12, margin = margin(r = 10)),
    axis.text.y = element_text(size = 10),
    plot.margin = margin(t=10, b=100, l = 10, r = 10)
  )

# Adjust variable contribution plot for PC2 with custom y-axis label
fviz_cos2_F_PC2 <- fviz_cos2(data_c1_F_scaled.pca, choice = "var", axes = 2, ggtheme = theme_bw(base_size = 10)) +
  ggtitle("Variable Contributions to PC2") +
  labs(y = "Contribution (%)") +  # Custom y-axis label
  scale_y_continuous(breaks = seq(0, 1, by = 0.1), labels = paste0(seq(0, 100, by = 10), "%")) +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.title.y = element_text(size = 12, margin = margin(r = 10)),
    axis.text.y = element_text(size = 10),
    plot.margin = margin(t=10, b=100, l = 10, r = 10)
  )

fviz_cos2_M_PC2 <- fviz_cos2(data_c1_M_scaled.pca, choice = "var", axes = 2, ggtheme = theme_bw(base_size = 10)) +
  ggtitle("Variable Contributions to PC2") +
  labs(y = "Contribution (%)") +  # Custom y-axis label
  scale_y_continuous(breaks = seq(0, 1, by = 0.1), labels = paste0(seq(0, 100, by = 10), "%")) +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.title.y = element_text(size = 12, margin = margin(r = 10)),
    axis.text.y = element_text(size = 10),
    plot.margin = margin(t=10, b=100, l = 10, r = 10)
  )




# Combine contribution plots for females
PCA_contribution_Female <- (fviz_cos2_F_PC1 | fviz_cos2_F_PC2) +
  plot_layout(heights = c(1, 1))

# Combine contribution plots for males
PCA_contribution_Male <- (fviz_cos2_M_PC1 | fviz_cos2_M_PC2) +
  plot_layout(heights = c(1, 1))


# Save plots
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/PCA_screeplot_Female.png", plot = PCA_screeplot_Female, width = 8, height = 6, units = "in", dpi = 300)

ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/PCA_contribution_Female.png", plot = PCA_contribution_Female, width = 8, height = 4, units = "in", dpi = 300)

ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/PCA_screeplot_Male.png", plot = PCA_screeplot_Male, width = 8, height = 6, units = "in", dpi = 300)

ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/PCA_contribution_Male.png", plot = PCA_contribution_Male, width = 8, height = 4, units = "in", dpi = 300)


#add PCs to the original data
val_c1_F_bin <- as.data.frame(data_c1_F_scaled.pca$x[,1:5])
data_c1_F_pca <- cbind(data_c1_F, val_c1_F_bin[1:5])

val_c1_M_bin <- as.data.frame(data_c1_M_scaled.pca$x[,1:5])
data_c1_M_pca <- cbind(data_c1_M, val_c1_M_bin[1:5])

Show code cell source Hide code cell source

# Convert images to base64
PCA_screeplot_Female <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/PCA_screeplot_Female.png")

PCA_contribution_Female <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/PCA_contribution_Female.png")

PCA_screeplot_Male <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/PCA_screeplot_Male.png")

PCA_contribution_Male <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/PCA_contribution_Male.png")

# Create the HTML
html_PCA_contribution_plots <- paste0("
<style>
  body, html {
    margin: 0; 
    padding: 0;
    /* If you want no horizontal scrollbar: */
    overflow-x: hidden; 
  }
  img {
    width: 100%; 
    height: auto; 
    display: block; 
    margin-bottom: 20px;
    border: 1px solid #ccc;
  }
</style>

<img src='", PCA_screeplot_Female, "' alt='Screeplot Female'>
<img src='", PCA_contribution_Female, "' alt='Contribution Plots Female'>
<img src='", PCA_screeplot_Male, "' alt='Screeplot Male'>
<img src='", PCA_contribution_Male, "' alt='Contribution Plots Male'>
")

# Display the HTML
IRdisplay::display_html(html_PCA_contribution_plots)

Show code cell source Hide code cell source

#visualize another way
PCA_Female <- fviz_pca_var(
  data_c1_F_scaled.pca, 
  col.var = "cos2",
  gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), 
  repel = TRUE
  ) +
  theme_bw(base_size = 10) +
  ggtitle("Females") +
  theme(
    axis.text.x = element_text(size = 12),
    axis.text.y = element_text(size = 12),
    axis.title.x = element_text(size = 14),
    axis.title.y = element_text(size = 14)
  )


PCA_Male <- fviz_pca_var(
  data_c1_M_scaled.pca, 
  col.var = "cos2",
  gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), 
  repel = TRUE
  ) +
  theme_bw(base_size = 10) +
  ggtitle("Males") +
  theme(
    axis.text.x = element_text(size = 12),
    axis.text.y = element_text(size = 12),
    axis.title.x = element_text(size = 14),
    axis.title.y = element_text(size = 14)
  )


# Save plots
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/PCA_Female.png", plot = PCA_Female, width = 8, height = 8, units = "in", dpi = 300)

ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/PCA_Male.png", plot = PCA_Male, width = 8, height = 8, units = "in", dpi = 300)

Show code cell source Hide code cell source

# Convert images to base64
PCA_Female <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/PCA_Female.png")

PCA_Male <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/PCA_Male.png")


# Create the HTML
html_PCA_plots <- paste0("
<style>
  body, html {
    margin: 0; 
    padding: 0;
    overflow-x: hidden; 
  }
  .image-row {
    display: flex;
    justify-content: space-between;
    align-items: center;
  }
  .image-row img {
    width: 49%; /* roughly 1/2 of the width; adjust if needed */
    height: auto;
    border: 1px solid #ccc;
  }
</style>
<div class='image-row'>
  <img src='", PCA_Female, "' alt='PCA Female'>
  <img src='", PCA_Male,   "' alt='PCA Male'>
</div>
")

IRdisplay::display_html(html_PCA_plots)

The first two principal components explain >85% of the data. These are the components that we will use for clustering.

---

4. K-means cluster analysis

K-means clustering is a machine learning method that groups data into a set number of clusters (k). It works by finding clusters where the points within each group are more similar to each other than to points in other groups.

The process begins by randomly placing cluster centers (centroids). Each data point is assigned to the closest centroid (usually based on Euclidean distance). Then, the centroids are updated to be the average position of all the points in their cluster. This process of assigning points and updating centroids repeats until the centroids stop moving significantly or a set number of steps is reached.

Cluster optimization

Finding the right number of clusters (k) is crucial in k-means clustering because it:

1. Improves model accuracy

   • Choosing an inappropriate k can lead to poorly defined clusters. Too few clusters may combine distinct groups, while too many can overfit the data.

2. Prevents Arbitrary Selection

   • Using systematic methods avoids guessing and ensures the chosen k is supported by the data.

3. Balances Complexity and Interpretability

   • A good k provides meaningful clusters without making the model overly complex or hard to interpret.

We use graphical methods to find the optimal k value in a k-means clustering algorithm. These methods ensure that k is chosen based on the structure of the data, leading to better and more reliable clustering outcomes. In this study, we will use the elbow method and silhouette method. The x-value at the inflection points indicate the optimal number of clusters.

Show code cell source Hide code cell source

# Determine optimal k value for clustering first two principal components

# Elbow method
elbow_Female <- fviz_nbclust(
  data_c1_F_pca[,11:12], kmeans, method = "wss") +
  geom_vline(xintercept = 2, linetype = 2) +
  theme_bw(base_size = 8) +
  ggtitle("A") +
  theme(
    plot.title = element_text(hjust = 0, size = 10),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    axis.title.x = element_text(size = 8),
    axis.title.y = element_text(size = 8)
  )
elbow_Male <- fviz_nbclust(
  data_c1_M_pca[,11:12], kmeans, method = "wss") +
  geom_vline(xintercept = 2, linetype = 2) +
  theme_bw(base_size = 8) +
  ggtitle("B") +
  theme(
    plot.title = element_text(hjust = 0, size = 10),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    axis.title.x = element_text(size = 8),
    axis.title.y = element_text(size = 8)
  )

# Silhouette method
silhouette_Female <- fviz_nbclust(
  data_c1_F_pca[,11:12], kmeans, method = "silhouette") +
  geom_vline(xintercept = 2, linetype = 2) +
  theme_bw(base_size = 8) +
  ggtitle("C") +
  theme(
    plot.title = element_text(hjust = 0, size = 10),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    axis.title.x = element_text(size = 8),
    axis.title.y = element_text(size = 8)
  )
silhouette_Male <- fviz_nbclust(
  data_c1_M_pca[,11:12], kmeans, method = "silhouette") +
  geom_vline(xintercept = 2, linetype = 2) +
  theme_bw(base_size = 8) +
  ggtitle("D") +
  theme(
    plot.title = element_text(hjust = 0, size = 10),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    axis.title.x = element_text(size = 8),
    axis.title.y = element_text(size = 8)
  )


# Combine plots for females and males with equal sizing
elbow_plots <- (elbow_Female | elbow_Male) +
  plot_layout(heights = c(1, 1)) +
  plot_annotation(
    theme = theme(plot.title = element_text(hjust = 0.5))
  )

silhouette_plots <- (silhouette_Female | silhouette_Male) +
  plot_layout(heights = c(1, 1)) +
  plot_annotation(
    theme = theme(plot.title = element_text(hjust = 0.5))
  )


# Save plots
ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/elbow_plots.png", plot = elbow_plots, width = 6, height = 3, units = "in", dpi = 300)

ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/silhouette_plots.png", plot = silhouette_plots, width = 6, height = 3, units = "in", dpi = 300)

Show code cell source Hide code cell source

# Convert images to base64
elbow_plots <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/elbow_plots.png")

silhouette_plots <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/silhouette_plots.png")



# Create the HTML
html_elbow_silhouette_plots <- paste0("
<style>
  body, html {
    margin: 0; 
    padding: 0;
    /* If you want no horizontal scrollbar: */
    overflow-x: hidden; 
  }
  img {
    width: 100%; 
    height: auto; 
    display: block; 
    margin-bottom: 20px;
    border: 1px solid #ccc;
  }
</style>

<img src='", elbow_plots, "' alt='elbow plots'>
<img src='", silhouette_plots, "' alt='silhouette plots'>
")

# Display the HTML
IRdisplay::display_html(html_elbow_silhouette_plots)

The graphs indicate that 2 to 4 clusters looks okay to try for females and males We can always change the number of clusters if it doesn’t make sense with our data later. Remember, clustering is a heuristic method (i.e., trial and error)!

Cluster analysis

k=2

Show code cell content Hide code cell content

#cluster data and add to dataset
km2_c1_F = kmeans(data_c1_F_pca[,11:12], centers = 2, iter.max = 100, nstart = 25)
km2_c1_F_bin <- cbind(data_c1_F_pca, Cluster=km2_c1_F$cluster)

km2_c1_M = kmeans(data_c1_M_pca[,11:12], centers = 2, iter.max = 100, nstart = 25)
km2_c1_M_bin <- cbind(data_c1_M_pca, Cluster=km2_c1_M$cluster)

Show code cell content Hide code cell content

# Plot clusters and PCA results
cb_palette_F <- c("indianred", "lightblue4")

labels <- c("bold(e[max])", "bold(Filter[max])", "bold(e[prop])", "bold(R[c])", "bold(G[c])", "bold(B[c])")

Female_km2 <- fviz_pca_biplot(
  data_c1_F_scaled.pca,
  habillage = as.factor(km2_c1_F_bin$Cluster),
  addEllipses = FALSE,
  mean.point = FALSE,
  palette = cb_palette_F,
  label = "none",
  pointshape = 16,
  pointsize = 2,
  repel = TRUE,
  col.var = "black"
) +
  labs(color = "Cluster") +
  theme_bw(base_size = 8) +
  ggtitle("A") +
  theme(
    plot.title = element_text(hjust = 0, size = 10),
    legend.position = "none",
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    axis.title.x = element_text(size = 8),
    axis.title.y = element_text(size = 8),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )
# Remove the default text layer if needed
Female_km2$layers <- Female_km2$layers[-2]

  # Add custom labels without lines
Female_km2_final <- Female_km2 + geom_text_repel(data = as.data.frame(get_pca_var(data_c1_F_scaled.pca)$coord),
                  aes(x = Dim.1, y = Dim.2, label = labels),
                  parse = TRUE,
                  size = 3,
                  segment.color = NA,
                  fontface = "bold")


cb_palette_M <- c("lightblue4", "indianred")

Male_km2 <- fviz_pca_biplot(
  data_c1_M_scaled.pca,
  habillage = as.factor(km2_c1_M_bin$Cluster),
  addEllipses = FALSE,
  mean.point = FALSE,
  palette = cb_palette_M,
  label = "none",
  pointshape = 16,
  pointsize = 2,
  repel = TRUE,
  col.var = "black"
) +
  labs(color = "Cluster") +
  theme_bw(base_size = 8) +
  ggtitle("B") +
  theme(
    plot.title = element_text(hjust = 0, size = 10),
    legend.position = "none",
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    axis.title.x = element_text(size = 8),
    axis.title.y = element_text(size = 8),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )
# Remove the default text layer if needed
Male_km2$layers <- Male_km2$layers[-2]

Male_km2_final <- Male_km2 + geom_text_repel(data = as.data.frame(get_pca_var(data_c1_M_scaled.pca)$coord),
                  aes(x = Dim.1, y = Dim.2, label = labels),
                  parse = TRUE,
                  size = 3,
                  segment.color = NA)


# Combine plots for females and males with equal sizing
km2_plots <- (Female_km2_final | Male_km2_final) +
  plot_layout(heights = c(1, 1)) +
  plot_annotation(
    theme = theme(plot.title = element_text(hjust = 0.5))
  )

ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/km2_plots.png", plot = km2_plots, width = 6, height = 3, units = "in", dpi = 300)

Show code cell source Hide code cell source

# Convert images to base64
km2_plots <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/km2_plots.png")



# Create the HTML
html_kmeans2_plots <- paste0("
<style>
  body, html {
    margin: 0; 
    padding: 0;
    /* If you want no horizontal scrollbar: */
    overflow-x: hidden; 
  }
  img {
    max-width: 800px;   /* ~8 inches at 100 dpi screen rendering */
    width: 100%;
    height: auto;
    display: block;
    margin-bottom: 20px;
    border: 1px solid #ccc;
  }
</style>

<img src='", km2_plots, "' alt='kmeans2 Plot'>
")

# Display the HTML
IRdisplay::display_html(html_kmeans2_plots)

Show code cell source Hide code cell source

#Fancier way of viewing clusters

#FEMALES 

# Define custom colors for clusters
custom_colors <- c("1" = "indianred", "2" = "lightblue4")

# Plot for females
pf_km2 <- plot_ly(km2_c1_F_bin, x = ~PC1, y = ~PC2, mode = "markers") %>%
  add_markers(
    size = 3.5,
    text = ~paste("BCPD:", BCPD, "<br>Morph_original:", Morph_original),
    color = ~as.character(Cluster),
    colors = custom_colors
  ) %>%
  layout(
    title = list(text = "Females (k=2)", font = list(size = 18)),
    xaxis = list(
      title = "PC1",
      showgrid = FALSE,
      zeroline = FALSE,
      showline = TRUE,
      linecolor = "black"
    ),
    yaxis = list(
      title = "PC2",
      showgrid = FALSE,
      zeroline = FALSE,
      showline = TRUE,
      linecolor = "black"
    ),
    plot_bgcolor = "white",  # White plot background like theme_bw()
    paper_bgcolor = "white", # White overall background
    font = list(family = "Arial", size = 10, color = "black")
  )


pf_km2

Show code cell content Hide code cell content

#Fancier way of viewing clusters

# MALES

# Define custom colors for clusters
custom_colors <- c("2" = "indianred", "1" = "lightblue4")

# Plot for males
pm_km2 <- plot_ly(km2_c1_M_bin, x = ~PC1, y = ~PC2, mode = "markers") %>%
  add_markers(
    size = 3.5,
    text = ~paste("BCPD:", BCPD, "<br>Morph_original:", Morph_original),
    color = ~as.character(Cluster),
    colors = custom_colors
  ) %>%
  layout(
    title = list(text = "Males (k=2)", font = list(size = 18)),
    xaxis = list(
      title = "PC1",
      showgrid = FALSE,
      zeroline = FALSE,
      showline = TRUE,
      linecolor = "black"
    ),
    yaxis = list(
      title = "PC2",
      showgrid = FALSE,
      zeroline = FALSE,
      showline = TRUE,
      linecolor = "black"
    ),
    plot_bgcolor = "white",  # White plot background like theme_bw()
    paper_bgcolor = "white", # White overall background
    font = list(family = "Arial", size = 10, color = "black")
  )


pm_km2

k=3

Show code cell content Hide code cell content

#cluster data and add to dataset
km3_c1_F = kmeans(data_c1_F_pca[,11:12], centers = 3, iter.max = 100, nstart = 25)
km3_c1_F_bin <- cbind(data_c1_F_pca, Cluster=km3_c1_F$cluster)

km3_c1_M = kmeans(data_c1_M_pca[,11:12], centers = 3, iter.max = 100, nstart = 25)
km3_c1_M_bin <- cbind(data_c1_M_pca, Cluster=km3_c1_M$cluster)

Show code cell content Hide code cell content

# Plot clusters and PCA results
cb_palette_F <- c("lightblue4", "indianred", "#117733")

labels <- c("bold(e[max])", "bold(Filter[max])", "bold(e[prop])", "bold(R[c])", "bold(G[c])", "bold(B[c])")

Female_km3 <- fviz_pca_biplot(
  data_c1_F_scaled.pca,
  habillage = as.factor(km3_c1_F_bin$Cluster),
  addEllipses = FALSE,
  mean.point = FALSE,
  palette = cb_palette_F,
  label = "none",
  pointshape = 16,
  pointsize = 2,
  repel = TRUE,
  col.var = "black"
) +
  labs(color = "Cluster") +
  theme_bw(base_size = 8) +
  ggtitle("A") +
  theme(
    plot.title = element_text(hjust = 0, size = 10),
    legend.position = "none",
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    axis.title.x = element_text(size = 8),
    axis.title.y = element_text(size = 8),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )
# Remove the default text layer if needed
Female_km3$layers <- Female_km3$layers[-2]

  # Add custom labels without lines
Female_km3_final <- Female_km3 + geom_text_repel(data = as.data.frame(get_pca_var(data_c1_F_scaled.pca)$coord),
                  aes(x = Dim.1, y = Dim.2, label = labels),
                  parse = TRUE,
                  size = 3,
                  segment.color = NA,
                  fontface = "bold")



cb_palette_M <- c("#117733", "indianred", "lightblue4")

Male_km3 <- fviz_pca_biplot(
  data_c1_M_scaled.pca,
  habillage = as.factor(km3_c1_M_bin$Cluster),
  addEllipses = FALSE,
  mean.point = FALSE,
  palette = cb_palette_M,
  label = "none",
  pointshape = 16,
  pointsize = 2,
  repel = TRUE,
  col.var = "black"
) +
  labs(color = "Cluster") +
  theme_bw(base_size = 8) +
  ggtitle("B") +
  theme(
    plot.title = element_text(hjust = 0, size = 10),
    legend.position = "none",
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    axis.title.x = element_text(size = 8),
    axis.title.y = element_text(size = 8),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )
# Remove the default text layer if needed
Male_km3$layers <- Male_km3$layers[-2]

  # Add custom labels without lines
Male_km3_final <- Male_km3 + geom_text_repel(data = as.data.frame(get_pca_var(data_c1_M_scaled.pca)$coord),
                  aes(x = Dim.1, y = Dim.2, label = labels),
                  parse = TRUE,
                  size = 3,
                  segment.color = NA,
                  fontface = "bold")


# Combine plots for females and males with equal sizing
km3_plots <- (Female_km3_final | Male_km3_final) +
  plot_layout(heights = c(1, 1)) +
  plot_annotation(
    theme = theme(plot.title = element_text(hjust = 0.5))
  )

ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/km3_plots.png", plot = km3_plots, width = 6, height = 3, units = "in", dpi = 300)

Show code cell source Hide code cell source

# Convert images to base64
km3_plots <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/km3_plots.png")



# Create the HTML
html_kmeans3_plots <- paste0("
<style>
  body, html {
    margin: 0; 
    padding: 0;
    /* If you want no horizontal scrollbar: */
    overflow-x: hidden; 
  }
  img {
    max-width: 800px;   /* ~8 inches at 100 dpi screen rendering */
    width: 100%;
    height: auto;
    display: block;
    margin-bottom: 20px;
    border: 1px solid #ccc;
  }
</style>

<img src='", km3_plots, "' alt='kmeans2 Plot'>
")

# Display the HTML
IRdisplay::display_html(html_kmeans3_plots)

Show code cell source Hide code cell source

#Fancier way of viewing clusters

# FEMALES

# Define custom colors for clusters
custom_colors <- c("1" = "lightblue4", "3" = "#117733", "2" = "indianred")

# Plot for females
pf_km3 <- plot_ly(km3_c1_F_bin, x = ~PC1, y = ~PC2, mode = "markers") %>%
  add_markers(
    size = 3.5,
    text = ~paste("BCPD:", BCPD, "<br>Morph_original:", Morph_original),
    color = ~as.character(Cluster),
    colors = custom_colors
  ) %>%
  layout(
    title = list(text = "Females (k=3)", font = list(size = 18)),
    xaxis = list(
      title = "PC1",
      showgrid = FALSE,
      zeroline = FALSE,
      showline = TRUE,
      linecolor = "black"
    ),
    yaxis = list(
      title = "PC2",
      showgrid = FALSE,
      zeroline = FALSE,
      showline = TRUE,
      linecolor = "black"
    ),
    plot_bgcolor = "white",  # White plot background like theme_bw()
    paper_bgcolor = "white", # White overall background
    font = list(family = "Arial", size = 10, color = "black")
  )

pf_km3

Show code cell source Hide code cell source

#Fancier way of viewing clusters

# MALES

# Define custom colors for clusters
custom_colors <- c("1" = "#117733", "3" = "lightblue4", "2" = "indianred")

# Plot for males
pm_km3 <- plot_ly(km3_c1_M_bin, x = ~PC1, y = ~PC2, mode = "markers") %>%
  add_markers(
    size = 3.5,
    text = ~paste("BCPD:", BCPD, "<br>Morph_original:", Morph_original),
    color = ~as.character(Cluster),
    colors = custom_colors
  ) %>%
  layout(
    title = list(text = "Males (k=3)", font = list(size = 18)),
    xaxis = list(
      title = "PC1",
      showgrid = FALSE,
      zeroline = FALSE,
      showline = TRUE,
      linecolor = "black"
    ),
    yaxis = list(
      title = "PC2",
      showgrid = FALSE,
      zeroline = FALSE,
      showline = TRUE,
      linecolor = "black"
    ),
    plot_bgcolor = "white",  # White plot background like theme_bw()
    paper_bgcolor = "white", # White overall background
    font = list(family = "Arial", size = 10, color = "black")
  )


pm_km3

k=4

Show code cell content Hide code cell content

#cluster data and add to dataset
km4_c1_F = kmeans(data_c1_F_pca[,11:12], centers = 4, iter.max = 100, nstart = 25)
km4_c1_F_bin <- cbind(data_c1_F_pca, Cluster=km4_c1_F$cluster)

km4_c1_M = kmeans(data_c1_M_pca[,11:12], centers = 4, iter.max = 100, nstart = 25)
km4_c1_M_bin <- cbind(data_c1_M_pca, Cluster=km4_c1_M$cluster)

Show code cell content Hide code cell content

# Plot clusters and PCA results
cb_palette_F <- c("darkorange", "#117733", "lightblue4", "indianred")

labels <- c("bold(e[max])", "bold(Filter[max])", "bold(e[prop])", "bold(R[c])", "bold(G[c])", "bold(B[c])")

Female_km4 <- fviz_pca_biplot(
  data_c1_F_scaled.pca,
  habillage = as.factor(km4_c1_F_bin$Cluster),
  addEllipses = FALSE,
  mean.point = FALSE,
  palette = cb_palette_F,
  label = "none",
  pointshape = 16,
  pointsize = 2,
  repel = TRUE,
  col.var = "black"
) +
  labs(color = "Cluster") +
  theme_bw(base_size = 8) +
  ggtitle("A") +
  theme(
    plot.title = element_text(hjust = 0, size = 10),
    legend.position = "none",
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    axis.title.x = element_text(size = 8),
    axis.title.y = element_text(size = 8),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )
# Remove the default text layer if needed
Female_km4$layers <- Female_km4$layers[-2]

  # Add custom labels without lines
Female_km4_final <- Female_km4 + geom_text_repel(data = as.data.frame(get_pca_var(data_c1_F_scaled.pca)$coord),
                  aes(x = Dim.1, y = Dim.2, label = labels),
                  parse = TRUE,
                  size = 3,
                  segment.color = NA,
                  fontface = "bold")

cb_palette_M <- c("indianred", "lightblue4", "#117733", "darkorange")

Male_km4 <- fviz_pca_biplot(
  data_c1_M_scaled.pca,
  habillage = as.factor(km4_c1_M_bin$Cluster),
  addEllipses = FALSE,
  mean.point = FALSE,
  palette = cb_palette_M,
  label = "none",
  pointshape = 16,
  pointsize = 2,
  repel = TRUE,
  col.var = "black"
) +
  labs(color = "Cluster") +
  theme_bw(base_size = 8) +
  ggtitle("B") +
  theme(
    plot.title = element_text(hjust = 0, size = 10),
    legend.position = "none",
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    axis.title.x = element_text(size = 8),
    axis.title.y = element_text(size = 8),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )
# Remove the default text layer if needed
Male_km4$layers <- Male_km4$layers[-2]

Male_km4_final <- Male_km4 + geom_text_repel(data = as.data.frame(get_pca_var(data_c1_M_scaled.pca)$coord),
                  aes(x = Dim.1, y = Dim.2, label = labels),
                  parse = TRUE,
                  size = 3,
                  segment.color = NA,
                  fontface = "bold")


# Combine plots for females and males with equal sizing
km4_plots <- (Female_km4_final | Male_km4_final) +
  plot_layout(heights = c(1, 1)) +
  plot_annotation(
    theme = theme(plot.title = element_text(hjust = 0.5))
  )

ggsave("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/km4_plots.png", plot = km4_plots, width = 6, height = 3, units = "in", dpi = 300)

Show code cell source Hide code cell source

# Convert images to base64
km4_plots <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/km4_plots.png")



# Create the HTML
html_kmeans4_plots <- paste0("
<style>
  body, html {
    margin: 0; 
    padding: 0;
    /* If you want no horizontal scrollbar: */
    overflow-x: hidden; 
  }
  img {
    max-width: 800px;   /* ~8 inches at 100 dpi screen rendering */
    width: 100%;
    height: auto;
    display: block;
    margin-bottom: 20px;
    border: 1px solid #ccc;
  }
</style>

<img src='", km4_plots, "' alt='kmeans2 Plot'>
")

# Display the HTML
IRdisplay::display_html(html_kmeans4_plots)

Show code cell source Hide code cell source

#Fancier way of viewing clusters

# FEMALES

# Define custom colors for clusters
custom_colors <- c("4" = "indianred", "3" = "lightblue4", "2" = "#117733", "1" = "darkorange")


# Plot for females
pf_km4 <- plot_ly(km4_c1_F_bin, x = ~PC1, y = ~PC2, mode = "markers") %>%
  add_markers(
    size = 3.5,
    text = ~paste("BCPD:", BCPD, "<br>Morph_original:", Morph_original),
    color = ~as.character(Cluster),
    colors = custom_colors
  ) %>%
  layout(
    title = list(text = "Females (k=4)", font = list(size = 18)),
    xaxis = list(
      title = "PC1",
      showgrid = FALSE,
      zeroline = FALSE,
      showline = TRUE,
      linecolor = "black"
    ),
    yaxis = list(
      title = "PC2",
      showgrid = FALSE,
      zeroline = FALSE,
      showline = TRUE,
      linecolor = "black"
    ),
    plot_bgcolor = "white",  # White plot background like theme_bw()
    paper_bgcolor = "white", # White overall background
    font = list(family = "Arial", size = 10, color = "black")
  )

pf_km4

#Fancier way of viewing clusters

# MALES

# Define custom colors for clusters
custom_colors <- c("4" = "darkorange", "3" = "#117733", "2" = "lightblue4", "1" = "indianred")

# Plot for males
pm_km4 <- plot_ly(km4_c1_M_bin, x = ~PC1, y = ~PC2, mode = "markers") %>%
  add_markers(
    size = 3.5,
    text = ~paste("BCPD:", BCPD, "<br>Morph_original:", Morph_original),
    color = ~as.character(Cluster),
    colors = custom_colors
  ) %>%
  layout(
    title = list(text = "Males (k=4)", font = list(size = 18)),
    xaxis = list(
      title = "PC1",
      showgrid = FALSE,
      zeroline = FALSE,
      showline = TRUE,
      linecolor = "black"
    ),
    yaxis = list(
      title = "PC2",
      showgrid = FALSE,
      zeroline = FALSE,
      showline = TRUE,
      linecolor = "black"
    ),
    plot_bgcolor = "white",  # White plot background like theme_bw()
    paper_bgcolor = "white", # White overall background
    font = list(family = "Arial", size = 10, color = "black")
  )


pm_km4

The boundaries between clusters are not well-defined, and our initial classification of color types does not consistently align with clear, distinct groups. As a result, evaluating the accuracy of our K-means clustering results is challenging without visual validation. To address this, we will export the clustering results and overlay them onto the outcomes of an additional analysis using additional R packages (see next section).

---

5. Export cluster analysis results

Add clusters to original dataset

# We need to add the clusters from clustering to the unstandardized dataset.
data_km2_c1_F_bin <- cbind(data_c1_F,Cluster=km2_c1_F_bin$Cluster)
data_km2_c1_M_bin <- cbind(data_c1_M,Cluster=km2_c1_M_bin$Cluster)

data_km3_c1_F_bin <- cbind(data_c1_F,Cluster=km3_c1_F_bin$Cluster)
data_km3_c1_M_bin <- cbind(data_c1_M,Cluster=km3_c1_M_bin$Cluster)

data_km4_c1_F_bin <- cbind(data_c1_F,Cluster=km4_c1_F_bin$Cluster)
data_km4_c1_M_bin <- cbind(data_c1_M,Cluster=km4_c1_M_bin$Cluster)


# To keep the dataset complete, we will add the morphs removed at the beginning (labelling them as "M")
#Label Morph_originals as cluster "M"
Cluster = c("M","M","M","M","M")
data_km2_c1_F_bin_morphs <- cbind(data_c1_F_morphs, Cluster)
Cluster = c("M","M","M","M","M","M","M","M","M","M","M","M","M","M","M")
data_km2_c1_M_bin_morphs <- cbind(data_c1_M_morphs, Cluster)

Cluster = c("M","M","M","M","M")
data_km3_c1_F_bin_morphs <- cbind(data_c1_F_morphs, Cluster)
Cluster = c("M","M","M","M","M","M","M","M","M","M","M","M","M","M","M")
data_km3_c1_M_bin_morphs <- cbind(data_c1_M_morphs, Cluster)

Cluster = c("M","M","M","M","M")
data_km4_c1_F_bin_morphs <- cbind(data_c1_F_morphs, Cluster)
Cluster = c("M","M","M","M","M","M","M","M","M","M","M","M","M","M","M")
data_km4_c1_M_bin_morphs <- cbind(data_c1_M_morphs, Cluster)


#Add Morphs to the dataset
data_km2_c1_F_final <- rbind(data_km2_c1_F_bin, data_km2_c1_F_bin_morphs)
data_km2_c1_M_final <- rbind(data_km2_c1_M_bin, data_km2_c1_M_bin_morphs)

data_km3_c1_F_final <- rbind(data_km3_c1_F_bin, data_km2_c1_F_bin_morphs)
data_km3_c1_M_final <- rbind(data_km3_c1_M_bin, data_km2_c1_M_bin_morphs)

data_km4_c1_F_final <- rbind(data_km4_c1_F_bin, data_km2_c1_F_bin_morphs)
data_km4_c1_M_final <- rbind(data_km4_c1_M_bin, data_km2_c1_M_bin_morphs)

Export data

# Export data and compare clusters to pictures to see if there are any obvious visual color patterns that fall out

write.csv(data_km2_c1_F_final,'C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/FEMALES/nomorphs_woutliers/rds_files/data_km2_c1_F_final.csv',row.names = TRUE)
write.csv(data_km2_c1_M_final,'C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/MALES/nomorphs_woutliers/rds_files/data_km2_c1_M_final.csv',row.names = TRUE)

write.csv(data_km3_c1_F_final,'C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/FEMALES/nomorphs_woutliers/rds_files/data_km3_c1_F_final.csv',row.names = TRUE)
write.csv(data_km3_c1_M_final,'C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/MALES/nomorphs_woutliers/rds_files/data_km3_c1_M_final.csv',row.names = TRUE)

write.csv(data_km4_c1_F_final,'C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/FEMALES/nomorphs_woutliers/rds_files/data_km4_c1_F_final.csv',row.names = TRUE)
write.csv(data_km4_c1_M_final,'C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/MALES/nomorphs_woutliers/rds_files/data_km4_c1_M_final.csv',row.names = TRUE)

---

6. Cluster evaluation

The R packages “recolorize” (Weller et al., 2024) and “patternize” (Belleghem et al., 2017) are tools that classify pixels into distinct color groups, which are then compared through PCA. We will use the results to (1) verify the optimal number of k-means clusters, and (2) assign individuals at the boundary between clusters to appropriate clusters.

Image alignment in patternize

We start with a folder of unaltered images of dorsal bodies. These images are exported PNGs of the measured dorsal body ROIs in ImageJ micatoolbox. Note that the images have slight variations in the size, shape, and angle, making it difficult to differentiate variation due to color pattern differences from that due to other factors. To address this, we set landmarks on the images to mark consistent reference points across all images, and perform alignment to standardize their orientation and scale. This ensures that subsequent analyses focus on meaningful color pattern differences rather than extraneous variability.

To perform image alignment in patternize, one set of XY coordinates of landmarks is required for each image. We do this in ImageJ using the multi-point tool and the custom “SetLandmarks.ijm” macro. Our landmarking scheme for the dorsal bodies has 9 points: 4 on the antipodal points of the body outline, and 1 on each midpoint of pereonites 2-6. The landmarks for each image should be exported as a two-column, tab-delimited text file with X coordinates on the left and Y coordinates on the right, and no header.

Once we have the landmark files, we can proceed with image alignment using the alignLan() function.

Females

# Set of specimen IDs
IDlist_c1_F <- tools::file_path_sans_ext(dir("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/FEMALES/nomorphs_woutliers/nobin/original_images/", ".png"))

# Make list with images
imageList_c1_F <- makeList(IDlist_c1_F, type = "image",
                      prepath = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/FEMALES/nomorphs_woutliers/nobin/original_images/",
                      extension = ".png")

# Make list with landmarks
landmarkList_c1_F <- makeList(IDlist_c1_F,
                         type = "landmark",
                         prepath = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/FEMALES/nomorphs_woutliers/nobin/landmarks/",
                         extension = "_landmarks.txt")

# Set target as BCPD_0102 (or whatever image you want)
target_c1_F <- landmarkList_c1_F[['BCPD_0102']]

# Set up mask
mask1_c1_F <- read.table("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/FEMALES/nomorphs_woutliers/nobin/masks/BCPD_0102_mask.txt", header = FALSE)

### Alignment ###
# This takes ~1 minute on a 16Gb RAM laptop running Ubuntu
imageList_aligned_c1_F <- alignLan(imageList_c1_F, landmarkList_c1_F, transformRef = target_c1_F,
                              adjustCoords = TRUE,
                              plotTransformed = FALSE,  # Suppress intermediate plots
                              resampleFactor = 5,
                              cartoonID = 'BCPD_0102',
                              maskOutline = mask1_c1_F)

# Save the aligned image list to an .rds file
saveRDS(imageList_aligned_c1_F, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/FEMALES/nomorphs_woutliers/rds_files/imageList_aligned_c1_F.rds")

Show code cell output Hide code cell output

[1] "sample 1 BCPD_0066 added to list"
[1] "sample 2 BCPD_0079 added to list"
[1] "sample 3 BCPD_0092 added to list"
[1] "sample 4 BCPD_0096 added to list"
[1] "sample 5 BCPD_0099 added to list"
[1] "sample 6 BCPD_0101 added to list"
[1] "sample 7 BCPD_0102 added to list"
[1] "sample 8 BCPD_0103 added to list"
[1] "sample 9 BCPD_0106 added to list"
[1] "sample 10 BCPD_0109 added to list"
[1] "sample 11 BCPD_0111 added to list"
[1] "sample 12 BCPD_0113 added to list"
[1] "sample 13 BCPD_0182 added to list"
[1] "sample 14 BCPD_0186 added to list"
[1] "sample 15 BCPD_0187 added to list"
[1] "sample 16 BCPD_0192 added to list"
[1] "sample 17 BCPD_0196 added to list"
[1] "sample 18 BCPD_0212 added to list"
[1] "sample 19 BCPD_0219 added to list"
[1] "sample 20 BCPD_0221 added to list"
[1] "sample 21 BCPD_0223 added to list"
[1] "sample 22 BCPD_0224 added to list"
[1] "sample 23 BCPD_0233 added to list"
[1] "sample 24 BCPD_0234 added to list"
[1] "sample 25 BCPD_0318 added to list"
[1] "sample 26 BCPD_0322 added to list"
[1] "sample 27 BCPD_0337 added to list"
[1] "sample 28 BCPD_0342 added to list"
[1] "sample 29 BCPD_0344 added to list"
[1] "sample 30 BCPD_0345 added to list"
[1] "sample 31 BCPD_0350 added to list"
[1] "sample 32 BCPD_0371 added to list"
[1] "sample 33 BCPD_0373 added to list"
[1] "sample 34 BCPD_0374 added to list"
[1] "sample 35 BCPD_0376 added to list"
[1] "sample 36 BCPD_0381 added to list"
[1] "sample 37 BCPD_0382 added to list"
[1] "sample 38 BCPD_0384 added to list"
[1] "sample 39 BCPD_0385 added to list"
[1] "sample 40 BCPD_0393 added to list"
[1] "sample 41 BCPD_0396 added to list"
[1] "sample 42 BCPD_0403 added to list"
[1] "sample 43 BCPD_0405 added to list"
[1] "sample 44 BCPD_0406 added to list"
[1] "sample 45 BCPD_0418 added to list"
[1] "sample 46 BCPD_0456 added to list"
[1] "sample 47 BCPD_0460 added to list"
[1] "sample 48 BCPD_0465 added to list"
[1] "sample 49 BCPD_0482 added to list"
[1] "sample 50 BCPD_0491 added to list"
[1] "sample 51 BCPD_0496 added to list"
[1] "sample 52 BCPD_0499 added to list"
[1] "sample 53 BCPD_0503 added to list"
[1] "sample 54 BCPD_0505 added to list"
[1] "sample 55 BCPD_0511 added to list"
[1] "sample 56 BCPD_0515 added to list"
[1] "sample 57 BCPD_0516 added to list"
[1] "sample 58 BCPD_0520 added to list"
[1] "sample 59 BCPD_0524 added to list"
[1] "sample 60 BCPD_0528 added to list"
[1] "sample 61 BCPD_0531 added to list"
[1] "sample 62 BCPD_0534 added to list"
[1] "sample 63 BCPD_0537 added to list"
[1] "sample 64 BCPD_0538 added to list"
[1] "sample 65 BCPD_0540 added to list"
[1] "sample 66 BCPD_0541 added to list"
[1] "sample 67 BCPD_0542 added to list"
[1] "sample 68 BCPD_0558 added to list"
[1] "sample 69 BCPD_0566 added to list"
[1] "sample 70 BCPD_0568 added to list"
[1] "sample 71 BCPD_0569 added to list"
[1] "sample 72 BCPD_0575 added to list"
[1] "sample 73 BCPD_0577 added to list"
[1] "sample 74 BCPD_0580 added to list"
[1] "sample 75 BCPD_0581 added to list"
[1] "sample 76 BCPD_0583 added to list"
[1] "sample 77 BCPD_0586 added to list"
[1] "sample 78 BCPD_0588 added to list"
[1] "sample 79 BCPD_0591 added to list"
[1] "sample 80 BCPD_0592 added to list"
[1] "sample 81 BCPD_0619 added to list"

[1] "sample 1 BCPD_0066 added to list"
[1] "sample 2 BCPD_0079 added to list"
[1] "sample 3 BCPD_0092 added to list"
[1] "sample 4 BCPD_0096 added to list"
[1] "sample 5 BCPD_0099 added to list"
[1] "sample 6 BCPD_0101 added to list"
[1] "sample 7 BCPD_0102 added to list"
[1] "sample 8 BCPD_0103 added to list"
[1] "sample 9 BCPD_0106 added to list"
[1] "sample 10 BCPD_0109 added to list"
[1] "sample 11 BCPD_0111 added to list"
[1] "sample 12 BCPD_0113 added to list"
[1] "sample 13 BCPD_0182 added to list"
[1] "sample 14 BCPD_0186 added to list"
[1] "sample 15 BCPD_0187 added to list"
[1] "sample 16 BCPD_0192 added to list"
[1] "sample 17 BCPD_0196 added to list"
[1] "sample 18 BCPD_0212 added to list"
[1] "sample 19 BCPD_0219 added to list"
[1] "sample 20 BCPD_0221 added to list"
[1] "sample 21 BCPD_0223 added to list"
[1] "sample 22 BCPD_0224 added to list"
[1] "sample 23 BCPD_0233 added to list"
[1] "sample 24 BCPD_0234 added to list"
[1] "sample 25 BCPD_0318 added to list"
[1] "sample 26 BCPD_0322 added to list"
[1] "sample 27 BCPD_0337 added to list"
[1] "sample 28 BCPD_0342 added to list"
[1] "sample 29 BCPD_0344 added to list"
[1] "sample 30 BCPD_0345 added to list"
[1] "sample 31 BCPD_0350 added to list"
[1] "sample 32 BCPD_0371 added to list"
[1] "sample 33 BCPD_0373 added to list"
[1] "sample 34 BCPD_0374 added to list"
[1] "sample 35 BCPD_0376 added to list"
[1] "sample 36 BCPD_0381 added to list"
[1] "sample 37 BCPD_0382 added to list"
[1] "sample 38 BCPD_0384 added to list"
[1] "sample 39 BCPD_0385 added to list"
[1] "sample 40 BCPD_0393 added to list"
[1] "sample 41 BCPD_0396 added to list"
[1] "sample 42 BCPD_0403 added to list"
[1] "sample 43 BCPD_0405 added to list"
[1] "sample 44 BCPD_0406 added to list"
[1] "sample 45 BCPD_0418 added to list"
[1] "sample 46 BCPD_0456 added to list"
[1] "sample 47 BCPD_0460 added to list"
[1] "sample 48 BCPD_0465 added to list"
[1] "sample 49 BCPD_0482 added to list"
[1] "sample 50 BCPD_0491 added to list"
[1] "sample 51 BCPD_0496 added to list"
[1] "sample 52 BCPD_0499 added to list"
[1] "sample 53 BCPD_0503 added to list"
[1] "sample 54 BCPD_0505 added to list"
[1] "sample 55 BCPD_0511 added to list"
[1] "sample 56 BCPD_0515 added to list"
[1] "sample 57 BCPD_0516 added to list"
[1] "sample 58 BCPD_0520 added to list"
[1] "sample 59 BCPD_0524 added to list"
[1] "sample 60 BCPD_0528 added to list"
[1] "sample 61 BCPD_0531 added to list"
[1] "sample 62 BCPD_0534 added to list"
[1] "sample 63 BCPD_0537 added to list"
[1] "sample 64 BCPD_0538 added to list"
[1] "sample 65 BCPD_0540 added to list"
[1] "sample 66 BCPD_0541 added to list"
[1] "sample 67 BCPD_0542 added to list"
[1] "sample 68 BCPD_0558 added to list"
[1] "sample 69 BCPD_0566 added to list"
[1] "sample 70 BCPD_0568 added to list"
[1] "sample 71 BCPD_0569 added to list"
[1] "sample 72 BCPD_0575 added to list"
[1] "sample 73 BCPD_0577 added to list"
[1] "sample 74 BCPD_0580 added to list"
[1] "sample 75 BCPD_0581 added to list"
[1] "sample 76 BCPD_0583 added to list"
[1] "sample 77 BCPD_0586 added to list"
[1] "sample 78 BCPD_0588 added to list"
[1] "sample 79 BCPD_0591 added to list"
[1] "sample 80 BCPD_0592 added to list"
[1] "sample 81 BCPD_0619 added to list"

[1] "sample 1 BCPD_0066 added to array"
[1] "sample 2 BCPD_0079 added to array"
[1] "sample 3 BCPD_0092 added to array"
[1] "sample 4 BCPD_0096 added to array"
[1] "sample 5 BCPD_0099 added to array"
[1] "sample 6 BCPD_0101 added to array"
[1] "sample 7 BCPD_0102 added to array"
[1] "sample 8 BCPD_0103 added to array"
[1] "sample 9 BCPD_0106 added to array"
[1] "sample 10 BCPD_0109 added to array"
[1] "sample 11 BCPD_0111 added to array"
[1] "sample 12 BCPD_0113 added to array"
[1] "sample 13 BCPD_0182 added to array"
[1] "sample 14 BCPD_0186 added to array"
[1] "sample 15 BCPD_0187 added to array"
[1] "sample 16 BCPD_0192 added to array"
[1] "sample 17 BCPD_0196 added to array"
[1] "sample 18 BCPD_0212 added to array"
[1] "sample 19 BCPD_0219 added to array"
[1] "sample 20 BCPD_0221 added to array"
[1] "sample 21 BCPD_0223 added to array"
[1] "sample 22 BCPD_0224 added to array"
[1] "sample 23 BCPD_0233 added to array"
[1] "sample 24 BCPD_0234 added to array"
[1] "sample 25 BCPD_0318 added to array"
[1] "sample 26 BCPD_0322 added to array"
[1] "sample 27 BCPD_0337 added to array"
[1] "sample 28 BCPD_0342 added to array"
[1] "sample 29 BCPD_0344 added to array"
[1] "sample 30 BCPD_0345 added to array"
[1] "sample 31 BCPD_0350 added to array"
[1] "sample 32 BCPD_0371 added to array"
[1] "sample 33 BCPD_0373 added to array"
[1] "sample 34 BCPD_0374 added to array"
[1] "sample 35 BCPD_0376 added to array"
[1] "sample 36 BCPD_0381 added to array"
[1] "sample 37 BCPD_0382 added to array"
[1] "sample 38 BCPD_0384 added to array"
[1] "sample 39 BCPD_0385 added to array"
[1] "sample 40 BCPD_0393 added to array"
[1] "sample 41 BCPD_0396 added to array"
[1] "sample 42 BCPD_0403 added to array"
[1] "sample 43 BCPD_0405 added to array"
[1] "sample 44 BCPD_0406 added to array"
[1] "sample 45 BCPD_0418 added to array"
[1] "sample 46 BCPD_0456 added to array"
[1] "sample 47 BCPD_0460 added to array"
[1] "sample 48 BCPD_0465 added to array"
[1] "sample 49 BCPD_0482 added to array"
[1] "sample 50 BCPD_0491 added to array"
[1] "sample 51 BCPD_0496 added to array"
[1] "sample 52 BCPD_0499 added to array"
[1] "sample 53 BCPD_0503 added to array"
[1] "sample 54 BCPD_0505 added to array"
[1] "sample 55 BCPD_0511 added to array"
[1] "sample 56 BCPD_0515 added to array"
[1] "sample 57 BCPD_0516 added to array"
[1] "sample 58 BCPD_0520 added to array"
[1] "sample 59 BCPD_0524 added to array"
[1] "sample 60 BCPD_0528 added to array"
[1] "sample 61 BCPD_0531 added to array"
[1] "sample 62 BCPD_0534 added to array"
[1] "sample 63 BCPD_0537 added to array"
[1] "sample 64 BCPD_0538 added to array"
[1] "sample 65 BCPD_0540 added to array"
[1] "sample 66 BCPD_0541 added to array"
[1] "sample 67 BCPD_0542 added to array"
[1] "sample 68 BCPD_0558 added to array"
[1] "sample 69 BCPD_0566 added to array"
[1] "sample 70 BCPD_0568 added to array"
[1] "sample 71 BCPD_0569 added to array"
[1] "sample 72 BCPD_0575 added to array"
[1] "sample 73 BCPD_0577 added to array"
[1] "sample 74 BCPD_0580 added to array"
[1] "sample 75 BCPD_0581 added to array"
[1] "sample 76 BCPD_0583 added to array"
[1] "sample 77 BCPD_0586 added to array"
[1] "sample 78 BCPD_0588 added to array"
[1] "sample 79 BCPD_0591 added to array"
[1] "sample 80 BCPD_0592 added to array"
[1] "sample 81 BCPD_0619 added to array"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0066 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0079 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0092 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0096 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0099 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0101 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0102 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0103 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0106 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0109 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0111 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0113 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0182 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0186 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0187 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0192 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0196 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0212 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0219 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0221 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0223 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0224 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0233 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0234 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0318 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0322 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0337 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0342 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0344 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0345 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0350 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0371 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0373 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0374 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0376 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0381 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0382 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0384 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0385 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0393 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0396 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0403 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0405 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0406 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0418 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0456 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0460 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0465 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0482 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0491 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0496 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0499 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0503 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0505 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0511 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0515 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0516 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0520 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0524 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0528 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0531 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0534 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0537 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0538 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0540 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0541 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0542 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0558 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0566 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0568 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0569 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0575 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0577 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0580 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0581 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0583 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0586 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0588 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0591 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0592 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0619 processed"

# Render combined plot
rasterstack_Female <- layout(matrix(1:50, nrow = 5)); par(mar = rep(1, 4))

# Set up the layout and margins
layout(matrix(1:50, nrow = 5))
par(mar = rep(1, 4))

# Plot each aligned image
lapply(imageList_aligned_c1_F, plotRasterstackAsImage)

Show code cell output Hide code cell output

$BCPD_0066

NULL

$BCPD_0079

NULL

$BCPD_0092

NULL

$BCPD_0096

NULL

$BCPD_0099

NULL

$BCPD_0101

NULL

$BCPD_0102

NULL

$BCPD_0103

NULL

$BCPD_0106

NULL

$BCPD_0109

NULL

$BCPD_0111

NULL

$BCPD_0113

NULL

$BCPD_0182

NULL

$BCPD_0186

NULL

$BCPD_0187

NULL

$BCPD_0192

NULL

$BCPD_0196

NULL

$BCPD_0212

NULL

$BCPD_0219

NULL

$BCPD_0221

NULL

$BCPD_0223

NULL

$BCPD_0224

NULL

$BCPD_0233

NULL

$BCPD_0234

NULL

$BCPD_0318

NULL

$BCPD_0322

NULL

$BCPD_0337

NULL

$BCPD_0342

NULL

$BCPD_0344

NULL

$BCPD_0345

NULL

$BCPD_0350

NULL

$BCPD_0371

NULL

$BCPD_0373

NULL

$BCPD_0374

NULL

$BCPD_0376

NULL

$BCPD_0381

NULL

$BCPD_0382

NULL

$BCPD_0384

NULL

$BCPD_0385

NULL

$BCPD_0393

NULL

$BCPD_0396

NULL

$BCPD_0403

NULL

$BCPD_0405

NULL

$BCPD_0406

NULL

$BCPD_0418

NULL

$BCPD_0456

NULL

$BCPD_0460

NULL

$BCPD_0465

NULL

$BCPD_0482

NULL

$BCPD_0491

NULL

$BCPD_0496

NULL

$BCPD_0499

NULL

$BCPD_0503

NULL

$BCPD_0505

NULL

$BCPD_0511

NULL

$BCPD_0515

NULL

$BCPD_0516

NULL

$BCPD_0520

NULL

$BCPD_0524

NULL

$BCPD_0528

NULL

$BCPD_0531

NULL

$BCPD_0534

NULL

$BCPD_0537

NULL

$BCPD_0538

NULL

$BCPD_0540

NULL

$BCPD_0541

NULL

$BCPD_0542

NULL

$BCPD_0558

NULL

$BCPD_0566

NULL

$BCPD_0568

NULL

$BCPD_0569

NULL

$BCPD_0575

NULL

$BCPD_0577

NULL

$BCPD_0580

NULL

$BCPD_0581

NULL

$BCPD_0583

NULL

$BCPD_0586

NULL

$BCPD_0588

NULL

$BCPD_0591

NULL

$BCPD_0592

NULL

$BCPD_0619

NULL

Males

# Set of specimen IDs
IDlist_c1_M <- tools::file_path_sans_ext(dir("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/MALES/nomorphs_woutliers/nobin/original_images/", ".png"))

# Make list with images
imageList_c1_M <- makeList(IDlist_c1_M, type = "image",
                      prepath = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/MALES/nomorphs_woutliers/nobin/original_images/",
                      extension = ".png")

# Make list with landmarks
landmarkList_c1_M <- makeList(IDlist_c1_M,
                         type = "landmark",
                         prepath = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/MALES/nomorphs_woutliers/nobin/landmarks/",
                         extension = "_landmarks.txt")

# Set target as BCPD_0392 (or whatever image you want)
target_c1_M <- landmarkList_c1_M[['BCPD_0392']]

# Set up mask
mask1_c1_M <- read.table("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/MALES/nomorphs_woutliers/nobin/masks/BCPD_0392_mask.txt", header = FALSE)

### Alignment ###
# This takes ~1 minute on a 16Gb RAM laptop running Ubuntu
imageList_aligned_c1_M <- alignLan(imageList_c1_M, landmarkList_c1_M, transformRef = target_c1_M,
                              adjustCoords = TRUE,
                              plotTransformed = FALSE,  # Suppress intermediate plots
                              resampleFactor = 5,
                              cartoonID = 'BCPD_0392',
                              maskOutline = mask1_c1_M)

# Save the aligned image list to an .rds file
saveRDS(imageList_aligned_c1_M, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/MALES/nomorphs_woutliers/rds_files/imageList_aligned_c1_M.rds")

Show code cell output Hide code cell output

[1] "sample 1 BCPD_0073 added to list"
[1] "sample 2 BCPD_0075 added to list"
[1] "sample 3 BCPD_0077 added to list"
[1] "sample 4 BCPD_0081 added to list"
[1] "sample 5 BCPD_0084 added to list"
[1] "sample 6 BCPD_0085 added to list"
[1] "sample 7 BCPD_0089 added to list"
[1] "sample 8 BCPD_0093 added to list"
[1] "sample 9 BCPD_0108 added to list"
[1] "sample 10 BCPD_0115 added to list"
[1] "sample 11 BCPD_0116 added to list"
[1] "sample 12 BCPD_0117 added to list"
[1] "sample 13 BCPD_0118 added to list"
[1] "sample 14 BCPD_0119 added to list"
[1] "sample 15 BCPD_0181 added to list"
[1] "sample 16 BCPD_0183 added to list"
[1] "sample 17 BCPD_0184 added to list"
[1] "sample 18 BCPD_0188 added to list"
[1] "sample 19 BCPD_0189 added to list"
[1] "sample 20 BCPD_0190 added to list"
[1] "sample 21 BCPD_0191 added to list"
[1] "sample 22 BCPD_0193 added to list"
[1] "sample 23 BCPD_0198 added to list"
[1] "sample 24 BCPD_0199 added to list"
[1] "sample 25 BCPD_0204 added to list"
[1] "sample 26 BCPD_0206 added to list"
[1] "sample 27 BCPD_0207 added to list"
[1] "sample 28 BCPD_0211 added to list"
[1] "sample 29 BCPD_0217 added to list"
[1] "sample 30 BCPD_0218 added to list"
[1] "sample 31 BCPD_0220 added to list"
[1] "sample 32 BCPD_0222 added to list"
[1] "sample 33 BCPD_0225 added to list"
[1] "sample 34 BCPD_0227 added to list"
[1] "sample 35 BCPD_0229 added to list"
[1] "sample 36 BCPD_0232 added to list"
[1] "sample 37 BCPD_0235 added to list"
[1] "sample 38 BCPD_0236 added to list"
[1] "sample 39 BCPD_0301 added to list"
[1] "sample 40 BCPD_0303 added to list"
[1] "sample 41 BCPD_0306 added to list"
[1] "sample 42 BCPD_0314 added to list"
[1] "sample 43 BCPD_0315 added to list"
[1] "sample 44 BCPD_0317 added to list"
[1] "sample 45 BCPD_0321 added to list"
[1] "sample 46 BCPD_0323 added to list"
[1] "sample 47 BCPD_0325 added to list"
[1] "sample 48 BCPD_0329 added to list"
[1] "sample 49 BCPD_0330 added to list"
[1] "sample 50 BCPD_0331 added to list"
[1] "sample 51 BCPD_0335 added to list"
[1] "sample 52 BCPD_0338 added to list"
[1] "sample 53 BCPD_0341 added to list"
[1] "sample 54 BCPD_0343 added to list"
[1] "sample 55 BCPD_0346 added to list"
[1] "sample 56 BCPD_0348 added to list"
[1] "sample 57 BCPD_0351 added to list"
[1] "sample 58 BCPD_0352 added to list"
[1] "sample 59 BCPD_0353 added to list"
[1] "sample 60 BCPD_0354 added to list"
[1] "sample 61 BCPD_0358 added to list"
[1] "sample 62 BCPD_0366 added to list"
[1] "sample 63 BCPD_0368 added to list"
[1] "sample 64 BCPD_0372 added to list"
[1] "sample 65 BCPD_0383 added to list"
[1] "sample 66 BCPD_0391 added to list"
[1] "sample 67 BCPD_0392 added to list"
[1] "sample 68 BCPD_0394 added to list"
[1] "sample 69 BCPD_0401 added to list"
[1] "sample 70 BCPD_0404 added to list"
[1] "sample 71 BCPD_0407 added to list"
[1] "sample 72 BCPD_0411 added to list"
[1] "sample 73 BCPD_0416 added to list"
[1] "sample 74 BCPD_0423 added to list"
[1] "sample 75 BCPD_0427 added to list"
[1] "sample 76 BCPD_0445 added to list"
[1] "sample 77 BCPD_0477 added to list"
[1] "sample 78 BCPD_0494 added to list"
[1] "sample 79 BCPD_0509 added to list"
[1] "sample 80 BCPD_0513 added to list"
[1] "sample 81 BCPD_0532 added to list"
[1] "sample 82 BCPD_0535 added to list"
[1] "sample 83 BCPD_0536 added to list"
[1] "sample 84 BCPD_0539 added to list"
[1] "sample 85 BCPD_0550 added to list"
[1] "sample 86 BCPD_0562 added to list"
[1] "sample 87 BCPD_0565 added to list"
[1] "sample 88 BCPD_0571 added to list"
[1] "sample 89 BCPD_0572 added to list"
[1] "sample 90 BCPD_0573 added to list"
[1] "sample 91 BCPD_0579 added to list"
[1] "sample 92 BCPD_0596 added to list"
[1] "sample 93 BCPD_0599 added to list"
[1] "sample 94 BCPD_0616 added to list"
[1] "sample 95 BCPD_0620 added to list"
[1] "sample 96 BCPD_0622 added to list"
[1] "sample 97 BCPD_0631 added to list"
[1] "sample 98 BCPD_0634 added to list"
[1] "sample 99 BCPD_0635 added to list"

[1] "sample 1 BCPD_0073 added to list"
[1] "sample 2 BCPD_0075 added to list"
[1] "sample 3 BCPD_0077 added to list"
[1] "sample 4 BCPD_0081 added to list"
[1] "sample 5 BCPD_0084 added to list"
[1] "sample 6 BCPD_0085 added to list"
[1] "sample 7 BCPD_0089 added to list"
[1] "sample 8 BCPD_0093 added to list"
[1] "sample 9 BCPD_0108 added to list"
[1] "sample 10 BCPD_0115 added to list"
[1] "sample 11 BCPD_0116 added to list"
[1] "sample 12 BCPD_0117 added to list"
[1] "sample 13 BCPD_0118 added to list"
[1] "sample 14 BCPD_0119 added to list"
[1] "sample 15 BCPD_0181 added to list"
[1] "sample 16 BCPD_0183 added to list"
[1] "sample 17 BCPD_0184 added to list"
[1] "sample 18 BCPD_0188 added to list"
[1] "sample 19 BCPD_0189 added to list"
[1] "sample 20 BCPD_0190 added to list"
[1] "sample 21 BCPD_0191 added to list"
[1] "sample 22 BCPD_0193 added to list"
[1] "sample 23 BCPD_0198 added to list"
[1] "sample 24 BCPD_0199 added to list"
[1] "sample 25 BCPD_0204 added to list"
[1] "sample 26 BCPD_0206 added to list"
[1] "sample 27 BCPD_0207 added to list"
[1] "sample 28 BCPD_0211 added to list"
[1] "sample 29 BCPD_0217 added to list"
[1] "sample 30 BCPD_0218 added to list"
[1] "sample 31 BCPD_0220 added to list"
[1] "sample 32 BCPD_0222 added to list"
[1] "sample 33 BCPD_0225 added to list"
[1] "sample 34 BCPD_0227 added to list"
[1] "sample 35 BCPD_0229 added to list"
[1] "sample 36 BCPD_0232 added to list"
[1] "sample 37 BCPD_0235 added to list"
[1] "sample 38 BCPD_0236 added to list"
[1] "sample 39 BCPD_0301 added to list"
[1] "sample 40 BCPD_0303 added to list"
[1] "sample 41 BCPD_0306 added to list"
[1] "sample 42 BCPD_0314 added to list"
[1] "sample 43 BCPD_0315 added to list"
[1] "sample 44 BCPD_0317 added to list"
[1] "sample 45 BCPD_0321 added to list"
[1] "sample 46 BCPD_0323 added to list"
[1] "sample 47 BCPD_0325 added to list"
[1] "sample 48 BCPD_0329 added to list"
[1] "sample 49 BCPD_0330 added to list"
[1] "sample 50 BCPD_0331 added to list"
[1] "sample 51 BCPD_0335 added to list"
[1] "sample 52 BCPD_0338 added to list"
[1] "sample 53 BCPD_0341 added to list"
[1] "sample 54 BCPD_0343 added to list"
[1] "sample 55 BCPD_0346 added to list"
[1] "sample 56 BCPD_0348 added to list"
[1] "sample 57 BCPD_0351 added to list"
[1] "sample 58 BCPD_0352 added to list"
[1] "sample 59 BCPD_0353 added to list"
[1] "sample 60 BCPD_0354 added to list"
[1] "sample 61 BCPD_0358 added to list"
[1] "sample 62 BCPD_0366 added to list"
[1] "sample 63 BCPD_0368 added to list"
[1] "sample 64 BCPD_0372 added to list"
[1] "sample 65 BCPD_0383 added to list"
[1] "sample 66 BCPD_0391 added to list"
[1] "sample 67 BCPD_0392 added to list"
[1] "sample 68 BCPD_0394 added to list"
[1] "sample 69 BCPD_0401 added to list"
[1] "sample 70 BCPD_0404 added to list"
[1] "sample 71 BCPD_0407 added to list"
[1] "sample 72 BCPD_0411 added to list"
[1] "sample 73 BCPD_0416 added to list"
[1] "sample 74 BCPD_0423 added to list"
[1] "sample 75 BCPD_0427 added to list"
[1] "sample 76 BCPD_0445 added to list"
[1] "sample 77 BCPD_0477 added to list"
[1] "sample 78 BCPD_0494 added to list"
[1] "sample 79 BCPD_0509 added to list"
[1] "sample 80 BCPD_0513 added to list"
[1] "sample 81 BCPD_0532 added to list"
[1] "sample 82 BCPD_0535 added to list"
[1] "sample 83 BCPD_0536 added to list"
[1] "sample 84 BCPD_0539 added to list"
[1] "sample 85 BCPD_0550 added to list"
[1] "sample 86 BCPD_0562 added to list"
[1] "sample 87 BCPD_0565 added to list"
[1] "sample 88 BCPD_0571 added to list"
[1] "sample 89 BCPD_0572 added to list"
[1] "sample 90 BCPD_0573 added to list"
[1] "sample 91 BCPD_0579 added to list"
[1] "sample 92 BCPD_0596 added to list"
[1] "sample 93 BCPD_0599 added to list"
[1] "sample 94 BCPD_0616 added to list"
[1] "sample 95 BCPD_0620 added to list"
[1] "sample 96 BCPD_0622 added to list"
[1] "sample 97 BCPD_0631 added to list"
[1] "sample 98 BCPD_0634 added to list"
[1] "sample 99 BCPD_0635 added to list"

[1] "sample 1 BCPD_0073 added to array"
[1] "sample 2 BCPD_0075 added to array"
[1] "sample 3 BCPD_0077 added to array"
[1] "sample 4 BCPD_0081 added to array"
[1] "sample 5 BCPD_0084 added to array"
[1] "sample 6 BCPD_0085 added to array"
[1] "sample 7 BCPD_0089 added to array"
[1] "sample 8 BCPD_0093 added to array"
[1] "sample 9 BCPD_0108 added to array"
[1] "sample 10 BCPD_0115 added to array"
[1] "sample 11 BCPD_0116 added to array"
[1] "sample 12 BCPD_0117 added to array"
[1] "sample 13 BCPD_0118 added to array"
[1] "sample 14 BCPD_0119 added to array"
[1] "sample 15 BCPD_0181 added to array"
[1] "sample 16 BCPD_0183 added to array"
[1] "sample 17 BCPD_0184 added to array"
[1] "sample 18 BCPD_0188 added to array"
[1] "sample 19 BCPD_0189 added to array"
[1] "sample 20 BCPD_0190 added to array"
[1] "sample 21 BCPD_0191 added to array"
[1] "sample 22 BCPD_0193 added to array"
[1] "sample 23 BCPD_0198 added to array"
[1] "sample 24 BCPD_0199 added to array"
[1] "sample 25 BCPD_0204 added to array"
[1] "sample 26 BCPD_0206 added to array"
[1] "sample 27 BCPD_0207 added to array"
[1] "sample 28 BCPD_0211 added to array"
[1] "sample 29 BCPD_0217 added to array"
[1] "sample 30 BCPD_0218 added to array"
[1] "sample 31 BCPD_0220 added to array"
[1] "sample 32 BCPD_0222 added to array"
[1] "sample 33 BCPD_0225 added to array"
[1] "sample 34 BCPD_0227 added to array"
[1] "sample 35 BCPD_0229 added to array"
[1] "sample 36 BCPD_0232 added to array"
[1] "sample 37 BCPD_0235 added to array"
[1] "sample 38 BCPD_0236 added to array"
[1] "sample 39 BCPD_0301 added to array"
[1] "sample 40 BCPD_0303 added to array"
[1] "sample 41 BCPD_0306 added to array"
[1] "sample 42 BCPD_0314 added to array"
[1] "sample 43 BCPD_0315 added to array"
[1] "sample 44 BCPD_0317 added to array"
[1] "sample 45 BCPD_0321 added to array"
[1] "sample 46 BCPD_0323 added to array"
[1] "sample 47 BCPD_0325 added to array"
[1] "sample 48 BCPD_0329 added to array"
[1] "sample 49 BCPD_0330 added to array"
[1] "sample 50 BCPD_0331 added to array"
[1] "sample 51 BCPD_0335 added to array"
[1] "sample 52 BCPD_0338 added to array"
[1] "sample 53 BCPD_0341 added to array"
[1] "sample 54 BCPD_0343 added to array"
[1] "sample 55 BCPD_0346 added to array"
[1] "sample 56 BCPD_0348 added to array"
[1] "sample 57 BCPD_0351 added to array"
[1] "sample 58 BCPD_0352 added to array"
[1] "sample 59 BCPD_0353 added to array"
[1] "sample 60 BCPD_0354 added to array"
[1] "sample 61 BCPD_0358 added to array"
[1] "sample 62 BCPD_0366 added to array"
[1] "sample 63 BCPD_0368 added to array"
[1] "sample 64 BCPD_0372 added to array"
[1] "sample 65 BCPD_0383 added to array"
[1] "sample 66 BCPD_0391 added to array"
[1] "sample 67 BCPD_0392 added to array"
[1] "sample 68 BCPD_0394 added to array"
[1] "sample 69 BCPD_0401 added to array"
[1] "sample 70 BCPD_0404 added to array"
[1] "sample 71 BCPD_0407 added to array"
[1] "sample 72 BCPD_0411 added to array"
[1] "sample 73 BCPD_0416 added to array"
[1] "sample 74 BCPD_0423 added to array"
[1] "sample 75 BCPD_0427 added to array"
[1] "sample 76 BCPD_0445 added to array"
[1] "sample 77 BCPD_0477 added to array"
[1] "sample 78 BCPD_0494 added to array"
[1] "sample 79 BCPD_0509 added to array"
[1] "sample 80 BCPD_0513 added to array"
[1] "sample 81 BCPD_0532 added to array"
[1] "sample 82 BCPD_0535 added to array"
[1] "sample 83 BCPD_0536 added to array"
[1] "sample 84 BCPD_0539 added to array"
[1] "sample 85 BCPD_0550 added to array"
[1] "sample 86 BCPD_0562 added to array"
[1] "sample 87 BCPD_0565 added to array"
[1] "sample 88 BCPD_0571 added to array"
[1] "sample 89 BCPD_0572 added to array"
[1] "sample 90 BCPD_0573 added to array"
[1] "sample 91 BCPD_0579 added to array"
[1] "sample 92 BCPD_0596 added to array"
[1] "sample 93 BCPD_0599 added to array"
[1] "sample 94 BCPD_0616 added to array"
[1] "sample 95 BCPD_0620 added to array"
[1] "sample 96 BCPD_0622 added to array"
[1] "sample 97 BCPD_0631 added to array"
[1] "sample 98 BCPD_0634 added to array"
[1] "sample 99 BCPD_0635 added to array"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0073 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0075 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0077 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0081 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0084 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0085 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0089 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0093 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0108 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0115 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0116 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0117 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0118 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0119 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0181 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0183 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0184 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0188 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0189 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0190 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0191 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0193 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0198 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0199 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0204 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0206 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0207 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0211 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0217 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0218 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0220 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0222 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0225 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0227 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0229 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0232 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0235 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0236 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0301 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0303 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0306 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0314 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0315 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0317 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0321 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0323 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0325 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0329 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0330 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0331 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0335 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0338 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0341 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0343 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0346 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0348 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0351 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0352 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0353 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0354 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0358 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0366 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0368 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0372 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0383 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0391 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0392 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0394 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0401 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0404 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0407 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0411 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0416 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0423 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0427 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0445 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0477 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0494 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0509 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0513 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0532 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0535 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0536 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0539 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0550 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0562 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0565 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0571 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0572 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0573 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0579 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0596 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0599 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0616 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0620 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0622 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0631 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0634 processed"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

Warning message:
"[rast] unknown extent"

[1] "sample BCPD_0635 processed"

# Render combined plot
rasterstack_Male <- layout(matrix(1:50, nrow = 5)); par(mar = rep(1, 4))


# Set up the layout and margins
layout(matrix(1:50, nrow = 5))
par(mar = rep(1, 4))

# Plot each aligned image
lapply(imageList_aligned_c1_M, plotRasterstackAsImage)

Show code cell output Hide code cell output

$BCPD_0073

NULL

$BCPD_0075

NULL

$BCPD_0077

NULL

$BCPD_0081

NULL

$BCPD_0084

NULL

$BCPD_0085

NULL

$BCPD_0089

NULL

$BCPD_0093

NULL

$BCPD_0108

NULL

$BCPD_0115

NULL

$BCPD_0116

NULL

$BCPD_0117

NULL

$BCPD_0118

NULL

$BCPD_0119

NULL

$BCPD_0181

NULL

$BCPD_0183

NULL

$BCPD_0184

NULL

$BCPD_0188

NULL

$BCPD_0189

NULL

$BCPD_0190

NULL

$BCPD_0191

NULL

$BCPD_0193

NULL

$BCPD_0198

NULL

$BCPD_0199

NULL

$BCPD_0204

NULL

$BCPD_0206

NULL

$BCPD_0207

NULL

$BCPD_0211

NULL

$BCPD_0217

NULL

$BCPD_0218

NULL

$BCPD_0220

NULL

$BCPD_0222

NULL

$BCPD_0225

NULL

$BCPD_0227

NULL

$BCPD_0229

NULL

$BCPD_0232

NULL

$BCPD_0235

NULL

$BCPD_0236

NULL

$BCPD_0301

NULL

$BCPD_0303

NULL

$BCPD_0306

NULL

$BCPD_0314

NULL

$BCPD_0315

NULL

$BCPD_0317

NULL

$BCPD_0321

NULL

$BCPD_0323

NULL

$BCPD_0325

NULL

$BCPD_0329

NULL

$BCPD_0330

NULL

$BCPD_0331

NULL

$BCPD_0335

NULL

$BCPD_0338

NULL

$BCPD_0341

NULL

$BCPD_0343

NULL

$BCPD_0346

NULL

$BCPD_0348

NULL

$BCPD_0351

NULL

$BCPD_0352

NULL

$BCPD_0353

NULL

$BCPD_0354

NULL

$BCPD_0358

NULL

$BCPD_0366

NULL

$BCPD_0368

NULL

$BCPD_0372

NULL

$BCPD_0383

NULL

$BCPD_0391

NULL

$BCPD_0392

NULL

$BCPD_0394

NULL

$BCPD_0401

NULL

$BCPD_0404

NULL

$BCPD_0407

NULL

$BCPD_0411

NULL

$BCPD_0416

NULL

$BCPD_0423

NULL

$BCPD_0427

NULL

$BCPD_0445

NULL

$BCPD_0477

NULL

$BCPD_0494

NULL

$BCPD_0509

NULL

$BCPD_0513

NULL

$BCPD_0532

NULL

$BCPD_0535

NULL

$BCPD_0536

NULL

$BCPD_0539

NULL

$BCPD_0550

NULL

$BCPD_0562

NULL

$BCPD_0565

NULL

$BCPD_0571

NULL

$BCPD_0572

NULL

$BCPD_0573

NULL

$BCPD_0579

NULL

$BCPD_0596

NULL

$BCPD_0599

NULL

$BCPD_0616

NULL

$BCPD_0620

NULL

$BCPD_0622

NULL

$BCPD_0631

NULL

$BCPD_0634

NULL

$BCPD_0635

NULL

Image segmentation in recolorize

To simplify the analysis, we reduce thousands of colors in the images to a manageable number for further refinement. This is achieved using the color histogram binning method in the recolorize package.

The color histogram binning method works by dividing each channel of a color space (e.g., RGB) into a predetermined number of bins. It then calculates the number of pixels falling into each bin and computes the average color for that bin. For instance, dividing each of the three RGB channels into two bins results in \(2^3 = 8\) total bins, with 8 representative colors.

Since patternize works with raster images and recolorize works with arrays, we first need to convert the raster objects into arrays to proceed with clustering.

After reducing the number of colors, we refine them by combining similar colors. This ensures the images are granular enough to capture important details while remaining simplified for analysis.

We use the recluster() function for this step, which:

1. Calculates the Euclidean distances between all color centers in a recolorize object.

2. Clusters them hierarchically using hclust.

3. Combines the most similar colors based on a user-defined cutoff.

To streamline the above steps, we use recolorize2(), a wrapper function that performs color histogram binning and refinement in one step. This function produces a list (rc_list) of simplified recolorize objects ready for further analysis.

By the end of this workflow, the images are reduced to a manageable set of representative colors while maintaining the necessary detail for meaningful comparisons.

# load the imageList_aligned object:
imageList_aligned_c1_F <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/FEMALES/nomorphs_woutliers/rds_files/imageList_aligned_c1_F.rds")

# convert from RasterBricks to image arrays:
imgs_F <- lapply(imageList_aligned_c1_F, brick_to_array)
names(imgs_F) <- names(imageList_aligned_c1_F)

# save raster extents:
extent_list_F <- lapply(imageList_aligned_c1_F, extent)


#remove former list, if needed. If we don't do this, you might get an error.
rm(c)


# fit initial recolorize fits (default color space is RGB). Note that recolorize2 is a wrapper function which runs recolorize and recluster sequentially in a single step.

rc_list_F <- lapply(imgs_F, 
                  function(i) recolorize2(i, bins = 3,
                                          cutoff = 45,
                                          plotting = FALSE,
                                          color_space = "Lab"))

# get all palettes and sizes
all_palettes_F <- do.call(rbind, lapply(rc_list_F, function(i) i$centers))
all_sizes_F <- do.call(c, lapply(rc_list_F, function(i) i$sizes / sum(i$sizes)))


# Next you can combine the color palettes from all of the recolorize objects in rc_list and use hclust_color to plot them and return a list of which colors to group together. 

# Note that as the number of clusters increases, the percent variation explained in the PCA decreases because your images are more complex (see chunk below for PCA).
cluster_list_F <- hclust_color(all_palettes_F, cutoff = 50, color_space="Lab")

# Save dendrogram to file
png(filename = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/dendrogram_Female.png", width = 800, height = 600, res = 150)

cluster_list_F <- hclust_color(all_palettes_F, cutoff = 50)

dev.off()

# make an empty matrix for storing the new palette
pod_palette_F <- matrix(NA, ncol = 3, nrow = length(cluster_list_F))

# for every color in cluster_list...
for (i in 1:length(cluster_list_F)) {
  
  # get the center indices
  idx <- cluster_list_F[[i]]
  
  # get the average value for each channel, using cluster size to get a weighted average
  ctr <- apply(all_palettes_F, 2, 
               function(j) weighted.mean(j[idx], 
                                         w = all_sizes_F[idx]))
  
  # store in the palette matrix
  pod_palette_F[i, ] <- ctr
}


# and apply
# set plotting = TRUE to view individual recolorize objects:
impose_list_F <- lapply(imgs_F, function(i) imposeColors(i, pod_palette_F, 
                                                     adjust_centers = FALSE, 
                                                     plotting = FALSE))


# save:
saveRDS(impose_list_F, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/FEMALES/nomorphs_woutliers/rds_files/recolorize_fits_c1_F.rds")


# convert back to patternize (including extent)
patternize_list_F <- lapply(impose_list_F, recolorize_to_patternize)
for (i in 1:length(patternize_list_F)) {
  for (j in 1:length(patternize_list_F[[1]])) {
    raster::extent(patternize_list_F[[i]][[j]]) <- extent_list_F[[i]]
  }
}

# and save
saveRDS(patternize_list_F, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/FEMALES/nomorphs_woutliers/rds_files/patternize_list_c1_F.rds")

Show code cell output Hide code cell output

Warning message in rm(c):
"object 'c' not found"

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

pdf: 2

# load the imageList_aligned object:
imageList_aligned_c1_M <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/MALES/nomorphs_woutliers/rds_files/imageList_aligned_c1_M.rds")

# convert from RasterBricks to image arrays:
imgs_M <- lapply(imageList_aligned_c1_M, brick_to_array)
names(imgs_M) <- names(imageList_aligned_c1_M)

# save raster extents:
extent_list_M <- lapply(imageList_aligned_c1_M, extent)



#remove former list, if needed. If we don't do this, you might get an error.
rm(c)


# fit initial recolorize fits (default color space is RGB). Note that recolorize2 is a wrapper function which runs recolorize and recluster sequentially in a single step.

rc_list_M <- lapply(imgs_M, 
                  function(i) recolorize2(i, bins = 3,
                                          cutoff = 45,
                                          plotting = FALSE,
                                          color_space = "Lab"))

# get all palettes and sizes
all_palettes_M <- do.call(rbind, lapply(rc_list_M, function(i) i$centers))
all_sizes_M <- do.call(c, lapply(rc_list_M, function(i) i$sizes / sum(i$sizes)))


# Next you can combine the color palettes from all of the recolorize objects in rc_list and use hclust_color to plot them and return a list of which colors to group together. 

# Note that as the number of clusters increases, the percent variation explained in the PCA decreases because your images are more complex (see chunk below for PCA).
cluster_list_M <- hclust_color(all_palettes_M, cutoff = 42)

# Save dendrogram to file
png(filename = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/dendrogram_Male.png", width = 800, height = 600, res = 150)

cluster_list_M <- hclust_color(all_palettes_M, cutoff = 42)

dev.off()

# make an empty matrix for storing the new palette
pod_palette_M <- matrix(NA, ncol = 3, nrow = length(cluster_list_M))

# for every color in cluster_list...
for (i in 1:length(cluster_list_M)) {
  
  # get the center indices
  idx <- cluster_list_M[[i]]
  
  # get the average value for each channel, using cluster size to get a weighted average
  ctr <- apply(all_palettes_M, 2, 
               function(j) weighted.mean(j[idx], 
                                         w = all_sizes_M[idx]))
  
  # store in the palette matrix
  pod_palette_M[i, ] <- ctr
}

# and apply
# set plotting = TRUE to view individual recolorize objects:
impose_list_M <- lapply(imgs_M, function(i) imposeColors(i, pod_palette_M, 
                                                     adjust_centers = FALSE, 
                                                     plotting = FALSE))

# save:
saveRDS(impose_list_M, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/MALES/nomorphs_woutliers/rds_files/recolorize_fits_c1_M.rds")


# convert back to patternize (including extent)
patternize_list_M <- lapply(impose_list_M, recolorize_to_patternize)
for (i in 1:length(patternize_list_M)) {
  for (j in 1:length(patternize_list_M[[1]])) {
    raster::extent(patternize_list_M[[i]][[j]]) <- extent_list_M[[i]]
  }
}

# and save
saveRDS(patternize_list_M, "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/MALES/nomorphs_woutliers/rds_files/patternize_list_c1_M.rds")

Show code cell output Hide code cell output

Warning message in rm(c):
"object 'c' not found"

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

Using 3^3 = 27 total bins

pdf: 2

Color pattern analyses in patternize

Since we now have the images segmented in the way that patternize needs, we can run any of the regular patternize functions on it (see the methods paper and examples repository). Here, we’ll use a custom function (https://hiweller.rbind.io/post/recolorize-patternize-workflow/) for running a PCA on the entire color pattern (all three colors simultaneously, rather than one color class at a time). This approach is similar to our initial PCA conducted before k-means clustering, but it uses binned colors rather than color pattern metrics derived from micaToolbox.

library(grid)

# load recolorize results
patternize_list_F <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/FEMALES/nomorphs_woutliers/rds_files/patternize_list_c1_F.rds")
patternize_list_M <- readRDS("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/MALES/nomorphs_woutliers/rds_files/patternize_list_c1_M.rds")


# load kmeans results
metadata_km2_F_c1 <- read.csv("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/FEMALES/nomorphs_woutliers/rds_files/data_km2_c1_F_final.csv")
metadata_km2_F_c1 <- metadata_km2_F_c1 %>%
  filter(Cluster != "M")
metadata_km2_M_c1 <- read.csv("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/MALES/nomorphs_woutliers/rds_files/data_km2_c1_M_final.csv")
metadata_km2_M_c1 <- metadata_km2_M_c1 %>%
  filter(Cluster != "M")

metadata_km3_F_c1 <- read.csv("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/FEMALES/nomorphs_woutliers/rds_files/data_km3_c1_F_final.csv")
metadata_km3_F_c1 <- metadata_km3_F_c1 %>%
  filter(Cluster != "M")
metadata_km3_M_c1 <- read.csv("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/MALES/nomorphs_woutliers/rds_files/data_km3_c1_M_final.csv")
metadata_km3_M_c1 <- metadata_km3_M_c1 %>%
  filter(Cluster != "M")

metadata_km4_F_c1 <- read.csv("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/FEMALES/nomorphs_woutliers/rds_files/data_km4_c1_F_final.csv")
metadata_km4_F_c1 <- metadata_km4_F_c1 %>%
  filter(Cluster != "M")
metadata_km4_M_c1 <- read.csv("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/MALES/nomorphs_woutliers/rds_files/data_km4_c1_M_final.csv")
metadata_km4_M_c1 <- metadata_km4_M_c1 %>%
  filter(Cluster != "M")


source("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/FEMALES/nomorphs_woutliers/rds_files/patPCA_total.R")

source("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/data/photos/PATTERNIZE/MALES/nomorphs_woutliers/rds_files/patPCA_total.R")


# run the PCA
pod_pca_F <- patPCA_total(patternize_list_F, quietly = FALSE)

pod_pca_M <- patPCA_total(patternize_list_M, quietly = FALSE)

# Reorder the rows of metadata files to match the order of the row names from pca files
metadata_km2_F_c1 <- metadata_km2_F_c1[match(rownames(pod_pca_F$x), metadata_km2_F_c1$BCPD), ]
metadata_km2_M_c1 <- metadata_km2_M_c1[match(rownames(pod_pca_M$x), metadata_km2_M_c1$BCPD), ]

metadata_km3_F_c1 <- metadata_km3_F_c1[match(rownames(pod_pca_F$x), metadata_km3_F_c1$BCPD), ]
metadata_km3_M_c1 <- metadata_km3_M_c1[match(rownames(pod_pca_M$x), metadata_km3_M_c1$BCPD), ]

metadata_km4_F_c1 <- metadata_km4_F_c1[match(rownames(pod_pca_F$x), metadata_km4_F_c1$BCPD), ]
metadata_km4_M_c1 <- metadata_km4_M_c1[match(rownames(pod_pca_M$x), metadata_km4_M_c1$BCPD), ]

# Calculate combined PCA limits for F and M
PCx <- 1; PCy <- 2
pca_limits <- list(
  x = range(c(pod_pca_F$x[, PCx], pod_pca_M$x[, PCx])),
  y = range(c(pod_pca_F$x[, PCy], pod_pca_M$x[, PCy]))
)

create_pca_plot <- function(pod_pca, rc_list, metadata, outline_colors, pca_limits, title) {
  # Set aspect ratio for slightly wider plots
  aspect_ratio <- 0.8  # Adjust to make the plot slightly wider

  pca_data <- as.data.frame(pod_pca$x[, c(PCx, PCy)])
  colnames(pca_data) <- c("PC1", "PC2")
  pca_data$image <- names(rc_list)
  pca_data$Cluster <- factor(metadata$Cluster)

  # Create base PCA plot
  plot <- ggplot(pca_data, aes(x = PC1, y = PC2)) +
    geom_point(size = 1, alpha = 0.8, aes(color = Cluster)) +
    scale_color_manual(values = outline_colors) +
    xlab(paste0("PC1 (", round(summary(pod_pca)$importance[2, PCx] * 100, 2), "% var.)")) +
    ylab(paste0("PC2 (", round(summary(pod_pca)$importance[2, PCy] * 100, 2), "% var.)")) +
    ggtitle(title) +
    theme_bw(base_size = 8) +
    theme(
      legend.position = "none",
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      axis.text.x = element_text(size = 8),
      axis.title.x = element_text(size = 8),
      axis.text.y = element_text(size = 8),
      axis.title.y = element_text(size = 8),
      plot.title = element_text(size = 10)
    ) +
    coord_fixed(ratio = aspect_ratio, xlim = pca_limits$x, ylim = pca_limits$y)  # Adjust aspect ratio

  # Adjust image size relative to PCA range
  image_size_x <- diff(pca_limits$x) * 0.05
  image_size_y <- image_size_x / aspect_ratio  # Adjust for aspect ratio

  for (i in 1:nrow(pca_data)) {
    img <- grid::rasterGrob(
      as.raster(recoloredImage(rc_list[[i]])),
      interpolate = TRUE
    )
    plot <- plot +
      annotation_custom(
        img,
        xmin = pca_data$PC1[i] - image_size_x / 2,
        xmax = pca_data$PC1[i] + image_size_x / 2,
        ymin = pca_data$PC2[i] - image_size_y / 2,
        ymax = pca_data$PC2[i] + image_size_y / 2
      )
  }

  # Adjust ellipse sizes to match image placement
  circle_radius_x <- image_size_x / 2
  circle_radius_y <- image_size_y / 2
  plot <- plot +
    ggforce::geom_ellipse(
      data = pca_data,
      aes(
        x0 = PC1, y0 = PC2,
        a = circle_radius_x, b = circle_radius_y, angle = 0, color = Cluster
      ),
      inherit.aes = FALSE,
      linewidth = 0.5  # Outline thickness
    )

  return(plot)
}

Show code cell output Hide code cell output

Summing raster lists...

Making dataframe from rasters...

Running PCA on 5 colors and 81 images...

done

Summing raster lists...

Making dataframe from rasters...

Running PCA on 5 colors and 99 images...

done

It can be hard to tell whether the PCA is capturing relevant axes of color pattern variation from a scatterplot, so we will instead use the actual images. We will also overlay the k-means clusters and see which cluster value (k=2, 3, or 4) shows the most reasonable clustering of morph types based on visual appearance.

Show code cell content Hide code cell content

# Define separate outline colors for each plot
outline_colors_F <- c("indianred4", "lightblue4", "#117733", "darkorange")  # Colors for females
names(outline_colors_F) <- c("1", "2", "3", "4")  # Adjust for your Cluster levels

outline_colors_M <- c("lightblue4", "indianred4", "#117733", "darkorange")  # Colors for males
names(outline_colors_M) <- c("1", "2", "3", "4")  # Adjust for your Cluster levels

# Create the PCA plots
recolorize_pca_Female_km2 <- create_pca_plot(pod_pca_F, rc_list_F, metadata_km2_F_c1, outline_colors_F, pca_limits, "C")
recolorize_pca_Female_km2 <- recolorize_pca_Female_km2

recolorize_pca_Male_km2 <- create_pca_plot(pod_pca_M, rc_list_M, metadata_km2_M_c1, outline_colors_M, pca_limits, "D")
recolorize_pca_Male_km2 <- recolorize_pca_Male_km2

recolorize_pca_Female_km3 <- create_pca_plot(pod_pca_F, rc_list_F, metadata_km3_F_c1, outline_colors_F, pca_limits, "C")
recolorize_pca_Female_km3 <- recolorize_pca_Female_km3

recolorize_pca_Male_km3 <- create_pca_plot(pod_pca_M, rc_list_M, metadata_km3_M_c1, outline_colors_M, pca_limits, "D")
recolorize_pca_Male_km3 <- recolorize_pca_Male_km3

recolorize_pca_Female_km4 <- create_pca_plot(pod_pca_F, rc_list_F, metadata_km4_F_c1, outline_colors_F, pca_limits, "C")
recolorize_pca_Female_km4 <- recolorize_pca_Female_km4

recolorize_pca_Male_km4 <- create_pca_plot(pod_pca_M, rc_list_M, metadata_km4_M_c1, outline_colors_M, pca_limits, "D")
recolorize_pca_Male_km4 <- recolorize_pca_Male_km4


# Combine plots for females and males with equal sizing
recolorize_pca_km2_plots <- (recolorize_pca_Female_km2 | recolorize_pca_Male_km2) +
  plot_layout(heights = c(1, 1)) +
  plot_annotation(
    theme = theme(plot.title = element_text(hjust = 0.5))
  )

recolorize_pca_km3_plots <- (recolorize_pca_Female_km3 | recolorize_pca_Male_km3) +
  plot_layout(heights = c(1, 1)) +
  plot_annotation(
    theme = theme(plot.title = element_text(hjust = 0.5))
  )

recolorize_pca_km4_plots <- (recolorize_pca_Female_km4 | recolorize_pca_Male_km4) +
  plot_layout(heights = c(1, 1)) +
  plot_annotation(
    theme = theme(plot.title = element_text(hjust = 0.5))
  )


# Save the plots with equal width and slightly adjusted height
ggsave(
  filename = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/recolorize_pca_km2_plots.png",
  plot = recolorize_pca_km2_plots,
  width = 6, height = 3, units = "in", dpi = 300
)

ggsave(
  filename = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/recolorize_pca_km3_plots.png",
  plot = recolorize_pca_km3_plots,
  width = 6, height = 3, units = "in", dpi = 300
)

ggsave(
  filename = "C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/recolorize_pca_km4_plots.png",
  plot = recolorize_pca_km4_plots,
  width = 6, height = 3, units = "in", dpi = 300
)

Show code cell source Hide code cell source

# Convert images to base64
recolorize_pca_km2_plots <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/recolorize_pca_km2_plots.png")

# Create the HTML
html_recolorize_pca_km2 <- paste0("
<style>
  body, html {
    margin: 0; 
    padding: 0;
    /* If you want no horizontal scrollbar: */
    overflow-x: hidden; 
  }
  img {
    max-width: 800px;   /* ~8 inches at 100 dpi screen rendering */
    width: 100%;
    height: auto;
    display: block;
    margin-bottom: 20px;
    border: 1px solid #ccc;
  }
</style>

<img src='", recolorize_pca_km2_plots, "' alt='recolorize kmeans2 Plot'>
")

IRdisplay::display_html(html_recolorize_pca_km2)


# Convert images to base64
recolorize_pca_km3_plots <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/recolorize_pca_km3_plots.png")

# Create the HTML
html_recolorize_pca_km3 <- paste0("
<style>
  body, html {
    margin: 0; 
    padding: 0;
    /* If you want no horizontal scrollbar: */
    overflow-x: hidden; 
  }
  img {
    max-width: 800px;   /* ~8 inches at 100 dpi screen rendering */
    width: 100%;
    height: auto;
    display: block;
    margin-bottom: 20px;
    border: 1px solid #ccc;
  }
</style>

<img src='", recolorize_pca_km3_plots, "' alt='recolorize kmeans3 Plot'>
")

IRdisplay::display_html(html_recolorize_pca_km3)


# Convert images to base64
recolorize_pca_km4_plots <- knitr::image_uri("C:/Users/bmc82/Documents/UF/PhD_Projects/DISSERTATION_MANUSCRIPT/Chapter_3/chapter3_data_analysis/images/recolorize_pca_km4_plots.png")

# Create the HTML
html_recolorize_pca_km4 <- paste0("
<style>
  body, html {
    margin: 0; 
    padding: 0;
    /* If you want no horizontal scrollbar: */
    overflow-x: hidden; 
  }
  img {
    max-width: 800px;   /* ~8 inches at 100 dpi screen rendering */
    width: 100%;
    height: auto;
    display: block;
    margin-bottom: 20px;
    border: 1px solid #ccc;
  }
</style>

<img src='", recolorize_pca_km4_plots, "' alt='recolorize kmeans4 Plot'>
")

IRdisplay::display_html(html_recolorize_pca_km4)

---

7. Summarize results

Camouflage types vary along a continuous spectrum of red, beige, and white. Cluster results plotted along the original principal components and the recolorize/patternize principal components show a discrete boundary between two continuous color types in males and females. Thus, Kmeans clustering with two clusters is the optimal clustering method for camouflage types. We will designate as red/beige (RB) and white/white-mottled (W).